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1.   INTRODUCTION 

For  the  man  on  the  street,  roundoff  analysis  of  floating- 
point computation  in  synonymous  with  keeping  bounds  on  accumulated 
errors:  at  any  stage  in  some  computation,  the  computed  approximation 
fl(A)  to  a  partial  result  A  satisfies 

fl(A)  =  A(l  +  6) 

where  6  =  (f 1 (A)  -  A)/A  is  the  relative  error  of  the  computation  [16]. 
The  idea  is  to  find  limits  on  the  magnitude  of  6. 

Recently,  however,  a  number  of  papers  have  appeared  adding 
probabilistic  considerations  to  the  distribution  of  roundoff  errors 
in  single  arithmetic  operations  and  of--what  almost  amounts  to  the  same- 
errors  in  representing  the  reals  as  floating-point  numbers  ([3],  [6], 
[14]).  They  show  that  relative  errors  are  not  uniformly  distributed 
within  their  bounds,  but  instead  cluster  around  zero.  The  error  bounds 
themselves  are  \/ery   pessimistic. 

In  this  paper  relative  errors  6  are  viewed  as  random  variables 
with  suitable  probability  distributions.  The  results  of  applying 
classical  probability  theory  to  roundoff  analysis  in  this  fashion  are: 
(1)  A  proof  of  the  generally  accepted  statement  that  accumulated  errors 
approach  normal  distributions  after  several  operations  ([13,  p.  113], 
[11,  p.  104],  [4,  p.  306]).  (2)  An  appreciation  of  the  rate  of  conver- 
gence to  normal  distributions  and  how  floating-point  addition  and  corre- 
lated errors  can  slow  this  convergence  down.  (3)  A  new  method  of 


automated  error  analysis  comparable  to  (but  in  most  ways  inferior  to) 
Interval  Analysis,  but  in  ewery   way  superior  to  Wilkinson  bound  analysis 

We  should  comment  that  it  is  not  "unrealistic"  to  treat  rela- 
tive errors  as  random  variables.  To  do  so  is  just  to  follow  the  basic 
philosophy  of  Bayesian  statistics:  if  you  don't  know  the  exact  value 
of  a  variable,  but  know  what  its  possible  values  are  and  what  the 
probability  of  taking  each  of  those  values  on  is,  then  you  should  treat 
the  variable  as  random  with  the  corresponding  probability  distribution. 
Treating  roundoff  errors  in  this  way  makes  our  work  fairly  clean,  plus 
lets  it  be  general  and  rigorous  as  well. 


2.  RELATIVE  ERRORS  AS  PROBABILITY  DENSITIES 

In  the  analysis  of  Wilkinson  [16]  the  relative  error 
6.  =  (fl(A)  -  A)/A  in  a  computed  result  A  is  rarely  known  exactly- 
usual  ly  one  says  that  it  lies  within  known  bounds.  These  bounds  are 
proportional  to  the  machine  precision  3  ,  where  3  is  the  arithmetic 
base  and  t  is  the  number  of  mantissa  digits  used  in  representing 
floating-point  numbers.  Computation  of  bounds  is  simple:  if  the 
floating-point  approximations  to  A  and  B  involve  the  relative  errors 
6.  and  6B,  respectively,  then  the  floating-point  approximation  to  A-B  is 

fl(A-B)  =  [A(l  +  6A)][B(1  +  6B)](1  +  6R) 

-  AB(1  +  6A  +  6B  +  6R  +  6A6B  +  8^R 

where  6D  is  the  rounding  error  incurred  in  the  multiplication  and  is 
bounded  by  the  machine  precision.  Since  hopefully  |6^|  «  1,  |6g|  «  1, 
we  have 

fl(AB)  -  AB(1  +  6A  +  6B  +  6R)  (2.2) 

to  a  yery   good  approximation.  Then  the  bound  for  the  error  on  fl(AB)  is 
the  sum  of  the  bounds  for  6«,  6g,  and  6R. 

Up  to  this  point  we  have  treated  each  relative  error  6  as  an 
unknown  element  of  a  known  interval.  It  is  not  a  big  step  to  view  6  as 
being  smeared  over  this  interval,  having  at  each  point  a  probability  of 
occurrence: 


Lower  Bound  0         Upper  Bound 

Thus  we  equate  6  with  an  underlying  probability  density  p(6)  so  that 


Probability  (5  £  x  ) 


°  p(x)  dx 


rX 


Lower  Bound 


p(x)  dx 


(2.3) 


Recall  that  a  probability  density  f  is  a  function  on  1R  satisfying 


(i)   f(x)  £  0    V  x  dR 


(ii)    f(x)  dx  =  1 

■  -00 


(2.4) 


and  that  the  probability  distribution  F  corresponding  to  the  density  f 
is  defined  by 

rX 


F(x) 


J-c 


f(t)  dt 


(2.5) 


The  question  arises  as  to  what  expressions  like  "A(l  +  6«)" 
mean  if  we  view  6«  as  a  probability  density.  This  is  resolved  by  the 
following  definition  and  explanation. 
Def :  A  smeared  number  (or  fuzzy  number)  is  a  pair  N  =  [A,f]  where  feF 

and  f  is  a  probability  density  on  1R.  The  probability  that  N  takes 

on  a  value  ^  A«C  for  any  CdR  is  given  by 


f  C 
Pr(N  i   A-C)  =   '  f(t)  dt 

J   -oo 


(2.6) 


The  idea  is  that  every  floating-point  approximation  to  a 
number  is  a  smeared  number.  Thus  fl(A)  =  A(l  +  6.)  =  [A,(l  +  OL 
We  clarify  this  concept  with  three  observations: 

Observation  1.  The  constant  1  in  expressions  of  the  form  (1  +  6) 

should  be  thought  of  as  an  impulse  function  A  (Dirac  delta) 
at  1. 


This  is  natural  since  1  should  represent  a  density  with  all  its 
support  at  1 .  In  the  following  we  will  use  1  and  A  interchangeably  when 
no  confusion  should  result. 

Observation  2.   If  6,  and  6~  are  independent  random  variables  with  corre- 
sponding densities  p,  and  Pp,  then  6  =  6,  +  6  corre- 
sponds to  the  convolution  density  p  =  p,  *  p?. 

Since  <5  is  the  density  which  gives  the  probability  at  each 
point  x  that  ^-i   +  <$?     =     x' 

Pr(6  *  x)     =  p^y)  p2(z)  dydz 

y  +  z  =  x 

•°°      fX    -    z 

p-|(y)  P2(z)  dydz 

J   -00      J   -CD 

Differentiating  with  respect  to  x, 


P(x)  = 


i-c 


p,  (x  -  z)  pJz)   dz  =  p,  *  p9(x) 


Note  from  observations  1  and  2  that  if  p.  is  the  density  corresponding 


to  6.,  then 


fl(A)  =  A(l  +  6A)  =  [A,  (1  +  6A)]  =  [A,  A*  pA] 

Observation  3.  If  6,  and  6~  are  independent  random  variables  with 
corresponding  densities  p,  and  p?,  then  6  =  <S,  6? 
corresponds  to  the  density  p  defined  by 
P(x)  = 


r°°  . 

Pl(x/z)  p2(z)  S|. 


The  derivation  is  similar  to  that  for  Observation  2.  We 
remark  that  the  approximation 

(1  +  6]  +  62  +  6^)  =  (1  +  6,  +  62) 

when  6-,  and  6?  are  small  random  variables  (with  corresponding  densities) 
still  makes  sense:  if  p,  is  zero  outside  [-D,  ,D..]  and  p„  is  zero  out- 
side [-D?,D?]  then  the  density  p  defined  as  above  for  6,60  is  zero  out- 
side [-D-|D?,D-,D?].  Assuming  D,  ,D?  «  1,  p  approximates  an  impulse 
function  at  zero  (the  identity  for  convolution). 

We  now  consider  what  fl (A  +  B)  looks  like.  Since 

fl(A  +  B)  =  (A(l  +  6A)  +  B(l  +  6A))  (1  +  6R) 

=  A(l  +  6A  +  6R)  +  B(l  +  6A  +  6R) 

where  again  6R  is  rounding  error,  we  start  with 

fl(A  +  B)  *  [A,  1  *  PA  *  PR]  +  [B,  1  *  pR  *  pR] 

The  results  for  subtraction  and  division  are  similar.  We  summarize  the 
results  here: 


fl 


(A  +  B)  =  [A,  1  *  pA  *  pR]  +  [B,  1  *  pB  *  pR]       (2.7) 


fl(A  -  B)  =  [A,  1  *  pA  *  pR]  +  [-B,  1  *  pB  *  pR]      (2.8) 

fl(A-B)    =  [A-B,  1  *  PA  *  PB  *  p  ]  (2.9) 

fl(A/B)    =  [A/B,  1  *  PA  *  pVB  *  pR]  (2.10) 

In  (2.10),  pvB  is  defined  by  pvg(x)  =  pfi(-x).  Notice  that  in  (2.7)  and 
(2.8)  the  result  is  not  actually  a  smeared  number  but  a  sum  of  smeared 
numbers.  This  reflects  Wilkinson's  statement  [16,  §28]  that  you  cannot 
say  much  about  the  error  in  a  sum  without  knowing  more  about  the  magni- 
tudes of  the  summands.  In  the  next  section  we  show  a  method  to  reduce 
the  sum  to  a  single  smeared  number.  In  short,  we  show  that  we  can  make 
a  consistent  algebra  for  smeared  numbers. 

We  must  now  select  probability  densities  p(6)  for  (1) 
representation  errors  (i.e.,  errors  incurred  by  representing  real 
numbers  in  floating-point  format)  and  (2)  the  rounding  error  6R  made 
in  single  arithmetic  operations.  All  other  errors  are  sums  of  these 
two  kinds.  As  mentioned  in  the  introduction,  we  can  use  the  same 
density  for  both  (1)  and  (2): 

Let  pR  be  the  uniform  probability  density  on 

[-61_t/2,   e1_t/2]  if  rounded  arithmetic  is  used  (2.11) 

[-6  "   ,  0]     if  chopped   (truncated)  arithmetic  is  used 


pD  for  Rounded  Arithmetic 

K 


pR  for  Chopped  Arithmetic 


Note  that  pR  is  a  uniform  distribution  inside  the  ordinary  Wilkinsonian 
bounds.  This  is  very  pessimistic  but  is  the  best  one  can  do  a  priori . 
If  we  consider  the  error  in  representing  A  =  a* 3  where  acp/3,1), 
e  e  H  then  since 

[fl(a)  -  a]  (2.12) 

-t    -t  -t 

is  uniformly  distributed  on  [-3  /2,3  /2]  for  rounding  or  [-3  ,0] 

for  chopping  (because  the  "distribution  of  (t  +  l)-th  digits"  is 

practically  uniform;  cf.  [14,  p.  270]) 


6A  =  (fl(A)  -  A)/A  =  (fl(a)  -  a)/a 


(2.13) 


is  uniformly  distributed  on  [-3"  /2a, 3~  /2a]  for  rounding  or  [-3"  /a,0] 


for  chopping.  If  we  assume  the  worst  case  a  =  1/3,  we  get  the  p 
defined  above. 


R 


3.  WHY  ACCUMULATED  ERROR  DENSITIES  ARE  ALMOST  NORMAL 

This  section  introduces  the  Central  Limit  Theorem,  which  is 
the  basis  for  the  popular  statement  that  accumulated  relative  errors 
become  normally  distributed.  The  Central  Limit  Theorem  is  useful 
intuitively  but  we  avoid  its  proof,  offering  Section  4  instead.  In 
this  section  we  explain  how  a  sum  of  relative  errors  converges  to  a 
normally  distributed  error.  After  showing  that  sums  of  smeared 
numbers  are  smeared  numbers  we  get  an  appreciation  of  how  this  con- 
vergence may  be  disturbed. 

We  begin  by  reviewing  some  important  concepts  from  probability 
theory.  Let  £  be  a  random  variable  having  corresponding  probability 
density,  /.  Then: 

Def:  The  nth  moment  M  of  E,   (or  of  /)  is 


f°° 


M. 


x*V(x)  dx.  (3.1) 


DefJ  Let  g  be  an  arbitrary  complex-valued  function  on  ]R.  Define  the 
expected  value  of  g(£)  by 


r°° 


E[g(5)]  = 


g(x)/(x)  dx  (3.2) 


Note  therefore  that  M   =  E[£n] 

n     L^  J 


Def:  The  mean  of  £  (or  of  /)  is  its  first  moment: 

y  =  E[£]  =  M]  (3.3) 
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Def :  The  variance  of  E,   (or  of  /)  is  given  by 
a2  =  E[(C  -  y)2]  =  M2  -  M2 

The  definitions  (3.1)  -  (3.4)  are  tied  together  with  the 
following  theorem: 


(3.4) 


Theorem  1 .  Given  densities  /, ,  ...,  /N  where  /.  has  mean  y.  and 

2 
variance  a.,  then  the  convolution  density  /,  *  ...  *  4     has  mean 

N    n  2    N   2       '        N 

y  =  I     y.  and  variance  a  s  5;  a.. 
i=l  n  i=l  n 

Proof.  Compute  the  Fourier  transforms  ("characteristic  functions," 

E[exp(iw£)])of  both  sides  of  /  =  /,  *  ...  *  /.,  and  use  the  facts  that 

oo   (_i)n  M 


(i)  g(w)  =  I 

n=0 


(i  =  /-I  ;  A  =  Fourier  Trans-  ) 
form 


(2)  g  *  h(w)  =  g(w)  h(w) 


Now  recall  that  the  Normal  (or  Gaussian)  probability  density  with 
parameters  u,a(a  >  0)  is 

tfx,  =  _J_e-(x-u)2/2o2 

/2tF  a 

so  its  probability  distribution  $(x)  is 


(3.5) 


*(x)  = 


1 


/2?  a  J 


rx   /.    x2/0  2 

e-(t  -  y)  /2a  dt 


(3.6) 


It  is  easy  to  show  that  the  mean  and  variance  of  \p   are  y  and  a  , 
respectively. 
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The  Central  Limit  Theorem  may  be  stated  as  follows:  given  a 

set  {6,,  ...,  6*.}   of  random  variables  having  respective  probability 

2 
densities  p.,  means  u.,  and  variances  a.  (i  =  1 ,  . . . ,  N)  then  under 

"certain  loose  conditions"  the  sum  variable  6  =  6-,  +  6?  +  . . .  +  6N 

(having  density  p  =  p,  *  p2  *  ...  *  pN)  converges  as  N  -*•  °°  to  a  normal 

2    2    2 
density  with  mean  u  =  u-.  +  vu  +  •  •  •  +  yN  and  variance  a  =  a-i  +  a2 

♦  ...  +  4 

Proof  of  the  Central  Limit  Theorem  is  dependent  of  course  on 
the  "loose  conditions"  surrounding  the  summand  variables  and  the  type 
of  convergence  used.  Good  surveys  may  be  found  in  [1]  and  [12].  We  will 
omit  any  proof  here  for  roundoff  errors  because  the  convergence  to  a 
normal  distribution,  although  useful  intuitively,  is  difficult  to  prove 
and  unnecessary  for  obtaining  statistical  bounds  (see  Section  4). 

The  idea  is  that  the  smeared  number  A(l  +6,  +  6?  +  ...  6N) 
is  equal  to  A(l  +  6)  where  6  has  a  density  that  is  almost  normal. 

Since: 

2 

1.  We  can  compute  fi's  mean  u  and  variance  a  easily. 

2.  Assuming  6  normal,  the  probability  that  |6  -  u|  =  3a 
is  0.002798  (i.e.,  a  99.7%  confidence  interval  on  6's 
value  is  [y  -  3a,  y  +  3a]). 

3.  For  N  larger  than  2  or  3  the  error  bound  [u  -  3a, 

u  +  3a]  is  inside  the  crude  Wilkinsonian  error  bound 
[-NB1_t/2,  N61-t/2]  for  Rounding  or  [-Ne1_t,  0]  for 
truncation 
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--apparently  without  much  effort  we  can  get  much  tighter,  more  realistic 
error  bounds  and  improve  our  error  estimates. 

This  gives  us  our  motivation:  statistical  error  bounds  seem 
tight.  But  we  need  more  details. 

We  compute  the  mean  and  variance  for  the  rounding  error  density 
pR  defined  by  (2.11).  For  rounding,  we  have 


^Rounding 


=  0 


"Rounding  =  ^"W/S  -  ^12  (3.7) 


For  chopping,  similarly 


Chopping    "B   /2 

a2  -     (e^)2  -   (-g1^)2     =     e2-2t  ,       . 

Chopping  3  2  p        /u  ^-i5) 


Therefore 


oD  dlf  aD    ,.    =  ar.    .    =  31-t//T2  (3.9) 

R       Rounding     Chopping       '  v   ' 

Notice  that  from  (2.9)  and  (2.10)  expressions  involving  only  repeated 

multipliation,  repeated  division,  or  mixed  multiplication  and  division 

5L___— — 

will  result  in  smeared  numbers  of  the  form  [A,  1  *  pR  *  pR  *  pR  *  ...  pR] 

where  the  number  m  of  pR's  in  the  convolution  product  depends  of  course 
on  the  number  of  operands  in  the  expression  and  the  accumulated  error 
for  each  operand.  For  example  if  A,B,C  are  irrational  then  fl(A(B*C)) 
■  [ABC,  1  *  pR  *  (pR  *  PR  *  PR)  *  PR]  since  each  number  has  a  representa- 
tion error  bounded  by  p  ,  and  the  two  multiplication  operations  have 

K 
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rounding  (truncation)  errors  bounded  by  pR.  The  Wilkinson  bounds  for  the 

5  1-t  5  1-t 
relative  error  in  ABC  would  be  6ABC  e  [--^  3   ,  j  3   ]  for  rounding  and 

^ABC  e  ^"^  "  '  ^  ^or  truncatl0n' 

If  we  assume  (Central  Limit  Theorem)  that  (pR  *  pR  *  pR  *  pR  *  pR) 

2    2  1-t 

is  normal  with  a  =  5aR,  y  =  0  for  rounding  or  y  =  5(-l/23   )  for  chopping, 

then  the  error  in  ABC  should  satisfy 

f  ["Vyf  ^1"t»Vrl  ^"^  Rounding 

^£[M-3o'y  +  3o]A(^f-3^-t,  (^s/^ 

Chopping 
99.7%  of  the  time.  Now  3^-rp-  *  1.94  <  5/2,  and  we  have  improved  our 
error  bounds  by  20%  after  only  two  operations.  As  the  number  of  operations 
N  increases  we  will  get  better  improvement  since  Wilkinson's  bound  are 
proportional  to  N  while  the  3a  bounds  are  proportional  to  •ft. 

It  is  not  unreasonable  to  assume  that  (pR  *  pR  *  pR  *  pR  *  pR) 
is  normal,  either.  Figures  1-4  show  the  result  of  repeated  convolution 
with  the  error  density  for  Rounding.  Already  (pR  *  pR  *  pR  *  pR)  is 
essentially  normal  (Figure  4).  On  a  superficial  basis  at  the  least, 
convergence  to  a  normal  density  seems  rapid. 

It  should  be  pointed  out  that  the  use  of  division  introduces 
"checked"  densities  pR  in  the  convolution  product  (c.f.  (2.10)).  This 
is  not  a  problem  because 

pR(x)  for  Rounding 

Pp(x)  =  Pp(-x)  ={  -,  t 

pR(x  -  6   )    for  Chopping      (3.10) 


14 


a   = 


61_t//T2" 


»t-l 


-26 


1-t 


0 


231_t    -231_t 


Figure  1 .  Figure  2. 

pR  for  Rounding  and  its  normal  approximation    Pd*Pr  and  normal  approximation 


o=B1-t/2 


■2B"  "  0  26 

Figure  3. 
pR*pR*pR  and  normal  approximation 


1-t 


a  =  ^'X//3 


Figure  4. 
p*p*p*pR  and  normal  approximation 
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Obviously  what  was  said  in  the  previous  paragraph  holds  without  modi- 
fication for  Rounding.  For  chopping  we  must  simply  be  more  careful  in 

v     1  1-t 
computation  of  the  mean  u,  since  the  mean  of  pR  is  +  «-3    instead  of 

-  iS        as  it  was  before.  For  Chopping,  division  errors  tend  to  cancel 

out  multiplication  errors  since  they  are  of  the  opposite  sign--this  is 

reflected  here  by  the  fact  that  division  makes  the  mean  of  the  final 

error  more  positive,  while  multiplication  makes  it  more  negative. 

We  commented  in  Section  2  that  addition  and  subtraction  were 
not  as  nice  as  multiplication  and  division  insofar  as  their  results 
were  not  easily  expressible  as  smeared  numbers.  For  the  rest  of  this 
section  we  address  the  problem  of  reducing  a  sum  of  smeared  numbers  to 
a  single  smeared  number. 

It  is  expedient  to  work  in  absolute  errors  instead  of  the 
relative  errors  we  have  been  working  with  until  now.  The  following 
definition  and  lemma  provide  the  necessary  groundwork. 

Def : For  a  bounded  density  /  and  nonzero  real  number  A,  the  stretched 
density  /  is  given  by 

/A(x)  ■  —/(f)  (3.1D 

|A|    A 

A 

We  say  that  /  is  /  stretched  by  the  value  A.  Here  we  have  restricted 
our  definition  to  bounded  densities  to  avoid  the  problem  introduced  by 
impulse  functions. 

Lemma  1  If  /  is  a  bounded  density,  then  [A,/]  =  [1,/  ]. 
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Proof  Since 


fA-C 


rA-C  1 


r(t)dt  =     ] 


^/(|)dt|y=t/A  =  (_  /(U)dy. 


[1,/  ]  is  equivalent  to  [A,/]  under  the  definition  of  smeared  numbers. 
Notice  that  [1,/  ]  is  an  "absolute"  smeared  number,  no  longer  relative 
to  A.  Effectively  [l,/A]  =  /A. 


Theorem  2  If  fl (A  +  B)  =  [A,  1  *  p  ]  +  [B,  1  *  p£]  and  A  f   0, 


B/0,  A  +  B  ,*  0  then  f  1(A  +  B)  =  [A  +  B,  1  *  p]  where  p(x) 
=  p{A/A+B)  *  p<B/A+B)(x) 


(3.12) 


Proof  Brute  force.  Note  that  (1  */)(x)  =/(x  -  1)  and  (1  */)  (x) 

*TK[flj -  *>■    Now: 

fl(A  +  B)     =     [A,    1   *  p,]   +  [B,   1    *  p9] 


=  [i,  (i  *  Pln  +  [i,  (i  *  p2n 
=  (i  *  p/  *.(i  *  p2)b 


(Lemma  1 


(This  follows  from  the  remark  that  [l,/]  =  /,  and  observation  2,  Part  2.) 

t=(A+B)y 


J   _oo  II  I   I 


A  +  B 


AB 


p2((A  *  B^  "  B)  % 


y=v+(B/A+B) 


A  +  B 
AB 


r°° 


,x  -  (A  +  B)v  -  (A  +  B) 

H*  A  ) 


/(A  +  B)vx  . 
P9(-i — 5 — H  dv 
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A  +  B 

L 

D(B/ 

A 

AB 
B 

A  +  B 
'A+B), 

A  +  B 

Pp               vv 

p(A/A+B)( 


A  +  B 


1   -  v) 


1 JA/A+B)  *  n(B/A+B),     x  n 

hTBTI  p2  lA  +  B  "    u 


A  +  B  H 


1         i   *  D(A/A+B)   *     (B/A+B),     x     , 
|A  +  B|  pl  p2  lA  +  BJ 

(1  *  p(A/A+B)  *  p(B/A+B)}A+B 


[1.  (1  *p|A/A+B)  *  p(B/  A+B))A+B] 


=  [A  +  B.  1  *  p{A/A+B)  *  4B/A+B)] 


Notice  that  A  -  B  =  A  +  (-B)  so  Theorem  2  applies  to  sub- 
traction also.  Theorem  2  leaves  us  a  few  steps  short  of  finishing  off 
all  the  problems  of  addition. 


Corollary  1 . 


n 


If       fl(_Z  A.)  =  [Ar  1  *  p1l  +  [A2,  1  *  p2]  + 

n         n 
then     fl(  J  A^  =  11     Ar  1  *  p] 


.  +  [A  ,  1  *  p  ] 
n     rn 


where 


1=1 


P  =  P 


i=l 
(A^ZA.)    (A0/^A.) 


T  i 


2'  i 


(A  /ZA.) 


n   i 


(3.13) 


(Proof.  Trivial.  Naturally  we  now  require  that  A.  f   0  (i  =  1 ,  . . . ,  n) 

n 

and  I  A.  f   0.) 
i=l  n 
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2 
Theorem  3.  Let  the  mean  and  variance  for  the  density  /  be  y  and  a  , 

a  2  2 

respectively.  Then  the  mean  and  variance  for  /  are  (Ay)  and  (A  a  ) , 

respectively.  (Proof.  Trivial.) 

From  coral lary  1  and  Theorem  3  it  is  immediately  obvious 

how  to  apply  the  central  limit  theorem  to  addition  and  subtraction. 

2 
If  the  relative  error  density  p,  of  A  has  mean  y,  and  variance  a1 

2 
and  the  density  p~  of  B  has  mean  y~  and  variance  a?,  then  the  density 

p  of  A  +  B  (as  in  Theorem  1)  has 

A  B 

mean        y  =  (A  +  B)y]  +  (A  +  B)  y2  (3.14) 

2     ,  A  \2  2   ,  B  ,2  2 


and  variance  a 


(/TTb)  al  +  (FT~B)  a2  (3J5) 


After  many  additions  the  error  density  will  converge  to  a  normal 

probability  density—but  more  operations  are  needed  because  the  stretched 

(A./ZA.) 
densities  p.  J   ""  will  not  be  identical  up  to  translation  as  they 

were  for  multiplication:  accumulated  errors  in  multiplication  become 
normally  distributed  much  more  quickly  than  those  in  addition.  The 
reader  who  has  researched  the  Central  Limit  Theorem  has  probably  dis- 
covered that  most  of  the  existing  proofs  assume  identically  distributed 
summand  variables.  Unidentical  distributions  can  slow  down  convergence-- 
especially  when  the  variables  differ  greatly.  Thus  when  stretching 
factors  like  A/(A  -  B)  became  large  (implying  lots  of  cancellation  and 
loss  of  significance  in  the  computation  of  A  -  B),  convergence  toward 
normal  error  densities  is  slowed. 

In  the  next  section  we  show  that  we  can  guarantee  confidence 
intervals  almost  as  good  as  [y  -  3a,  u  +  3a]  all  of  the  time,  even  with 
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slowed  convergence.  Thus  convergence  and  the  Central  Limit  Theorem 
in  general  should  be  used  only  for  an  intuitive  feel  of  what  is  going 
on,  while  the  growth  of  a  is  worth  serious  attention.  Note  that  when 
stretching  factors  like  A/(A  -  B)  get  large,  convergence  not  only  slows 
down  but  variances  get  large  also  (Theorem  3).  This  is  a  serious 
weakness  of  statistical  error  analysis. 
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4.  PROOF  THAT  STATISTICAL  BOUNDS  HOLD  FOR  ROUNDOFF  ERRORS 

Section  3  introduced  the  Central  Limit  Theorem,  showed  that  it 
could  be  used  to  show  how  accumulated  roundoff  errors  have  asymptotically 
normal  distributions,  and  discussed  the  rate  of  convergence  to  normal 
distributions.  The  literature  surrounding  the  Central  Limit  Theorem 
is  labyrinthine  since  there  are  any  number  of  assumptions  one  can  make 
about  the  summand  variables  which  simplify  or  aggravate  the  proof.  How- 
ever, as  we  mentioned  earlier  we  are  not  so  much  concerned  with  proving 
convergence  as  being  able  to  say,  for  any  roundoff  variable  6,  some- 
thing like 

Pr( |6  -  y|  *  3a)  ^  .0028  (4.1) 

as  we  can  with  normally  distributed  random  variables.  In  general  given 
a  fixed  small  probability  p  we  want  to  find  a  corresponding  X  such  that 

Pr(|5  -  u|  =  Xa)     -     P.  (4.2) 

Since  6  will  correspond  to  a  lengthy  convolution  of  error  densities, 
determining  X   exactly  will  be  unfeasible.  Instead  we  can  get  bounds 
on  A.  It  would  be  simpler  to  just  assume  that  6  is  normally  distributed 
and  look  up  A  in  normal  distribution  tables.  However,  this  is  not 
always  legitimate:  see  e.g.  [5].  Although  three  or  four  convolutions 
can  approximate  a  normal  density  well  (see  Figures  1-4)  large  stretching 
factors  can  upset  this  approximation.  The  purpose  of  this  section  is  to 
bound  how  much  these  factors  can  upset  it,  by  obtaining  bounds  for  A  in 
(4.2).  We  start  with  Chebyshev's  Inequality. 
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Theorem  4.  Let  6  be  a  random  variable  with  underlying  density  p(t), 

2 
having  mean  0  and  variance  a  .  Then 

Pr(  1 6 |  =  Aa)  =  1/A2 


Proof,  a2  =  E[62] 


J_c 


r  p(t)  dt 


=>  L> t2p(t)  dt 

P(t)  dt 


'|t|=Ac 

(Aa)2 


|t|=Aa 


=  (Aa)2  Pr(|6|  =  Aa) 


Pr ( | 6 1  =  Aa)  =  -j   ,  as  desired, 


A 


This  says  that  if  we  want 


Pr(|6|  =  Aa)  =  .0028 


then  we  will  be  guaranteed  this  relationship  provided 


1 


0028  =s>   A  =  18.89 


Clearly  this  is  not  very  good  since  A  =  3  works  for  normally  distributed 

variables,  and  after  many  convolutions  we  expect  6  to  be  very   close  to 

normal.  It  turns  out  we  can  obtain  a  much  tighter  limit  on  A  if  we  go 

through  the  analysis  used  to  obtain  the  Chernoff  bound. 

N 
Theorem  5.  Let  6  =  I   6.  be  a  roundoff  error,  i.e.  s<5  is  represented  by 

i=l  n 
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the  convolution  of  N  (possibly  stretched)  copies  of  pR,  the  rounding 
density.  Then 


Pr(|6|  =  Xo)     =     2e 


2 
where  a  is  the  variance  of  6. 


Proof.  Set  d(x)  =  {°   J  f  *° 
Then 


d(x)  =  e(x)  —  exp(L(x-Xa))  for  any  positive  L, 


This  implies 

Pr(6  =  \o)     =  E[d(6)] 
=  E[e(6)] 

=  E[exp(L6)  exp(-LAa)] 
=  E[exp(L6)]  e"LAa 

Now  if  p(t)  is  the  density  associated  with  6  then 


E[exp(L6)]     = 


exp(Lt)  p(t)  dt 

J   -00 

00  f]         fCO 


I     — j-    t  p(t)  dt  [p  is  of  compact  support] 
n=0  n-  J-" 


oo  m  Ln 

=  y  J] 

n=0 

where  M  is  the  nth  moment  of  p(t).  Since  p(t)  is  symmetric  about  0 
(i.e.,  even)  the  odd  terms  vanish  and  we  are  left  with 
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oo  M   L2n 
E[exp(LS)]  =  I       2" 
n=0  K      '' 

We  now  need  the  following  lemma. 

Lemma  2.  For  roundoff  densities  p(t) 

M2n   <  (2n)l 
(M2)n     n!2n 

Proof.  It  can  be  shown  that  if  p(t)  =  p,  *  p~  *  ...  *  PN(t)  and  M  . 

is  the  nth  moment  of  p.  ,  that  the  nth  moment  M  of  p  is 

k  n 

n        ,      N 

Mn  =  I   n  'n  '     n  '  , n  Mn  k 
n     Q  nrn2.  ...  nN.  k=]  nk,K 

n1+...+nN=n 


This  is  done  as  in  Theorem  1  by  comparing  the  characteristic  functions 

(essentially  Fourier  transforms)  of  both  sides  of  the  equation 

N 
p  =  II  *  p  .  We  then  use  the  fact  that  the  p,  are  stretched  versions 

k=l   k  k 

of  the  roundoff  density  pR  (2.11),  with  moments 


[Ak(e1_t/2)]n 

n,k 


n  even 


0  n  odd 

(A  ) 

where  p.  =  pD   ,  i.e.,  A.  is  a  nonzero  stretching  factor.  Pulling  all 
k    K  k 

of  this  material  together  leads  to  the  bound  of  the  lemma. 

It  is  interesting  that  if  p(t)  is  normal  with  mean  0, 
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(M2)n    a2"    n!  2n  (M2  =  °  ' 

--so  we  have  equality  with  the  stated  bound.  Therefore: 

1.  The  bound  is  fairly  tight,  especially  as  N  gets  large. 

2.  The  lemma  holds  for  other  densities  besides  stretched 
versions  of  pR. 

3.  There  is  probably  an  elegant  proof  of  this  lemma. 
It  now  follows  that 

oo     L  M  2  2 

E[exp(L6)]  =  I     ]-r   (--T^)n  =  exp(^-f-) 

n=0  n'    c  c 

and  therefore 

>     <     L2  2 
P(6  =  Xo)     =  exp( — y —  LAa)  for  all  positive  L. 

The  Chernoff  Bound  is  now  obtained  by  choosing  L  such  that  the  right 
hand  side  of  this  expression  is  optimal.   By  differentiating  it  is  easy 
to  see  that  L  =  A/a  is  the  correct  choice. 
It  now  follows  immediately  that 

2 

Pr(5  =  Xo)     =     e"A  /2  and 

2 

Pr(|6|  =  Xo)     =     2e"X  /2  as  stated. 

The  Chernoff  Bound  on  the  probability  is  much  tighter  than  the 
one  obtained  by  Chebyshev's  Inequality,  and  is  almost  as  good  as  having 
a  normal  distribution.  For  comparison  we  tabulate  the  A's  obtained  from 
each  method: 
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Chernoff  Bound  for    Chebyshev  Bound  for 
>  X  X  X 

Pr(|6|  =  Xo)  (6  is  normal)  (6  is  roundoff  error)     (6  arbitrary) 


2.448  3.162 

2.716  4.472 
3.035  7.072 
3.255  10.00 
3.625  18.90 

3.717  22.36 
3.899  31.62 
4.292  70.71 
4.799  223.6 
5.257  707.2 


10 

1.645 

05 

1.960 

02 

2.326 

01 

2.576 

002798 

3 

002 

3.090 

001 

3.290 

0002 

3.719 

00002 

4.265 

000002 

4.753 

Table  1.  Bounds  for  X  given  information  about  6 

Table  1  shows  us  that  99.7%  confidence  intervals  on  the  magnitude  of  6 
are  (-3.625a,  3.625a);  this  result  compares  quite  closely  with  the 
(-3a,  +  3a)  bounds  we  could  use  if  it  were  known  that  6  were  normal, 
and  has  been  used  in  the  program  that  actually  does  statistical  error 
analysis  (Section  6) . 
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5.  CORRELATED  ERRORS 

Until  now  we  have  been  assuming  that  all  our  relative  errors 
were  independent  random  variables.  Now  we  consider  the  case  in  which 
they  are  not  all  independent,  i.e.,  some  of  the  errors  are  related  to 
others  in  a  common  expression.  This  is  the  case,  for  example,  in  the 
evaluation  of  fl(A-A).  Assuming  that 


fl(A)  =  A(l  +  6A)  =  [A,l  *  pA]  (5.1) 

then  our  prior  analysis  would  claim 

fl(A-A)  =  A-A(l  +  5A  +  6A  +  6R) 

=  [A-A.1  *  pA  *  pA  *  pR].  (5.2) 

This  is  not  correct,  however;  it  is  too  optimistic.  The  correct  ex- 
pression is 


fl(A-A)  =  A-A(l  +  26A  +  6R) 


=  [A-A.A  *  p^2)  *  pR]  (5.3) 


(2) 
where  pi  '  is  pA  stretched  by  the  factor  2,  as  in  (3.11);  i.e., 

pA2)(t)  =  lPA(t/2)-  (5-4) 

We  remark  that  the  rounding  error  6R  (with  density  pR)  is  always  initially 
an  independent  random  variable.  Correlation  comes  from  repeated  use  of 
an  existing  error  as  in  the  case  with  6.  above;  new  rounding  errors  are 
always  uncorrelated  with  anything. 
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The  generalization  of  the  above  analysis  to  the  general  case 

is  straightforward.  Consider  a  smeared  number  A 

(c-,)    (c?)        (cN) 
A  =  [A,l  *  P]    *  p2    *  ...  *  PN  N  ]  (5.5) 

where  p.  is  a  probability  density  (k  =  1,  ...,  N)  and  the  C.  are 
stretching  factors.  Suppose  that  two  distinct  error  densities  p., 
p.  correspond  to  the  same  error  variable  6.  Without  loss  of  generality 
(convolution  is  commutative)  we  can  assume  these  densities  are   p,  and 


p?;  A  can  then  be  rewritten  as 

(Ci  +  < 
A  =  [A,l  *  p1  '      *  p3  J  *  ...  *  PN  n   ]         (5.6) 


(c,  +  c0)         (c0)        (cN) 


This  corresponds  to  merging  the  two  references  to  6  into  one.  Other 
references  to  6  should  also  be  merged. 

Note  that,  because  convolution  of  densities  (or  addition  of 
random  variables)  is  commutative,  correlated  errors  must  be  merged 
throughout  computation.  For  example,  in  computing 

X  =  fl(---(((A  +  B)  +  C)  +  D)  +  E)  +  ...  +Z)  +  A)      (5.7) 

the  error  6.  involved  in  the  last  step  of  the  addition  will  be  slightly 
correlated  with  the  error  of  the  rest  of  the  sum;  the  resulting  smeared 
number  is 

X  -  [X,  1  *  P<2A/X)  *  P<B/X)  *  ...  *P<Z/X)]         (5.8) 
where  X=2A+B+C+...+Z. 
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6.  THE  SNORT  SYSTEM  AND  PROBLEMS  OF  IMPLEMENTATION 

As  this  project  of  statistical  error  analysis  progressed,  it 
was  hoped  that  the  final  results  could  be  automated  to  perform  a_ 
posteriori  error  bound  calculation  in  the  spirit  of  Interval  Analysis. 
Since  for  small  relative  errors  6  =  (fl(x)  -  x)/x  we  have 

x  z     fl(x)  (1  -  6) 

it  seems  natural  to  evaluate  fl(x)  in  the  natural  way  and  use  statisti- 
cal bounds  on  6  to  bound  the  error  in  fl(x).  While  Interval  Analysis 
will  give  results  like 

x  e     [xlow'  xhigh] 
we  would  give  results  like 

x  e  [fl(x)  (1  +  u  -  B),  fl(x)  (1  +  u  +  B)] 

2 
where  B  =  (3.625a)  is  the  99.7%-confidence  bound  on  6,  u  and  a  being 

the  computed  mean  and  variance  of  6  viewed  as  a  random  variable. 

Automation  of  the  statistical  analysis  was  in  fact  completed  by 
development  of  SNORTRAN,  a  PL/I-like  language.  It  is  similar  in 
philosophy  to  PL/I-FORMAC  in  that  SNORTRAN  source  is  preprocessed  and 
expanced  into  PL/I,  which  is  then  handled  by  the  PL/I  compiler.  The 
preprocessor,  a  SNOBOL  routine  included  at  the  end  of  this  section,  is 
unsophisticated:  current  need  did  not  warrant  sophistication. 

The  results  of  several  test  programs  (also  included  at  the 
end  of  this  section)  reflect  negatively  on  continued  use  of  statistical 
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relative  error  analysis  for  almost  all  problems.  There  are  at  least 
four  reasons  for  the  discouragement: 

1.  There  is  no  way  to  represent  0  as  a  smeared  number. 

2.  In  subtraction  of  two  numbers  with  similar  values,  can- 
cellation causes  enormous  increase  in  error  variance 
(stretching  factors  get  large).  One  or  two  such  can- 
cellations will  result  in  a  value  of  a  so  large  that 
"useless"  error  bounds  result. 

3.  In  only  the  most  well-conditioned,  stable  computations 
(e.g.,  repeated  multiplication)  are  statistical  bounds 
competitive  with  Interval  Analysis. 

4.  Correlation  between  errors  cannot  be  easily  taken  into 
account  (as  shown  in  Section  5)  so  statistical  "bounds" 
ignoring  correlation  may  use  smaller  stretching  factors 
than  they  should,  and  give  incorrectly  small  results. 

The  actual  implementation  uses  Interval  Analysis  methods  to 
generate  worst-case  values  of  a  and  y,  hence  statistical  bounds  gen- 
erated by  the  program  should  be  bounds  (ignoring  correlation  effects) 
There  are  four  basic  SNORT  commands,  all  of  them  keyword-driven  for 

simplicity: 

TMTTTni  T7r  j  ROUNDED) 
INITIALIZE  {cHOPPEDJ' 

SMEAR  <variable>  [,<variable>]; 
EVAL  <variable>  =  <expression>; 
PRINT  <variable>  [,<variable>]; 
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INITIALIZE  causes  the  preprocessor  to  copy  SNORT  system 
declarations  and  routines  into  the  translated  program  output,  depending 
on  the  type  of  arithmetic  being  simulated.  If  rounded  arithmetic  is 
being  used,  the  preprocessor  copies  in  the  file  SNORT. RND,  otherwise 
copying  SNORT. CHP.  The  former  is  reproduced  at  the  end  of  this  section; 
the  latter,  requiring  full  Interval  Analysis  to  compute  error  means,  was 
never  written. 

SMEAR  <variable>  results  in  declarations  for  <variable>  which 
permit  it  to  be  used  as  a  smeared  number.  "<variable>"  may  be  sub- 
scripted if  desired  (e.g.,  SMEAR  A  (10,0:5))  and  a  list  of  <variable>s 
generates  a  list  of  appropriate  declarations. 

EVAL  <variable>  =  <expression>  causes  an  SLR(1 )  parser  to 
break  down  <expression>  into  SNORT  system  subroutine  calls.  Currently 
<expression>  should  only  involve  the  operators  +,-,/,*  although  this 
could  easily  be  expanded  or  modified  by  changing  the  parser  tables.  All 
variables  in  <variable>  and  <expression>  should  be  SMEARed--the  pre- 
processor does  not  currently  check. 

PRINT  provides  formatted  output  of  smeared  variables.  It 
lists  computed  values,  3.625a  bounds  and  the  value  bounds  they  imply, 
and  Interval  Analysis  value  bounds. 

Four  sample  runs  are  included  in  the  following  pages,  in- 
volving algorithms  which  highlight  progressively  more  serious  problems 
with  statistical  relative  error  analysis.  The  first  shows  that  even 
in  the  best  case  of  serial  multiplication,  statistical  analysis  does 
not  give  better  bounds  than  interval  analysis.  The  other  three  show 
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respectively  what  happens  when  addition,  subtraction,  and  computed  zeroes 
are  added. 

Preliminary  example:     Computation  of  sin  (1.0) 

The  first  test  of  the  SNORTRAN  system  attempted  to  compute 
sin(l)  =  .84147098...  via  the  little-used  relationship 


9      n=l     2n 


This  concentrates  on  the  roundoff  properties  of  multiplication 
(ignoring  computation  of  the  cosines)  which  we  know  statistical  analysis 
handles  well.  Unfortunately  the  results  are  not  encouraging.  Statisti- 
cal analysis  obtains  the  bounds 

sin(l)  e  (.84146956,  .84147235) 

after  32  factors,  while  Interval  Analysis  gets 

sin(l)  e  (.8414694,  .8414710). 


32 


» 
< 

LU 

CD 


** 

«    H 

* 

* 

* 

#    ♦ 

». 

t 

1 

Z 

t 

1 

1    •• 

• 

• 

•     1 

—   CC 

— • 

Of 

•     1 

#  0 

•o 

c 

•     1 

»  »- 

-4 

QC 

# 

If.     | 

0  0 

V. 

or 

CO 

vO     | 

Ul   < 

X 

LLI 

a: 

«n    1 

O  LL 

4~ 

□ 

0  1 

O       1 

Cl 

z 

t- 

in    1 

O  CO 

a. 

X 

a 

O 

1 

O  O 

1-4 

u 

— 

< 

-0    I 

O  O 

I 

i- 

LL 

0>    | 

O 

CO 

— 

< 

00    | 

O  # 

z 

CD 

t- 

LU 

i-    1 

O 

c 

V 

2: 

z 

O    1 

O  1- 

mi 

x 

UJ 

1-4 

1 

O  *- 

t- 

■~- 

1/1 

l/"> 

CO    1 

O  03 

< 

i/, 

LU 

0 

-*  1 

^ 

O       | 

_j 

c 

or 

0 

CO     I 

at 

O  < 

LU 

o 

a. 

<*  1 

0  0 

0T 

UJ 

(M 

0  1 

H 

O  a. 

*  » 

— 

a 

i-O 

1 

uoo 

or 

LU 

sT 

r-   1 

< 

O     I 

Lli 

X 

•^ 

0 

a. 

^   ! 

9W 

a.  0  z 

a 

K 

X 

z 

LU 

-  1 

tt 

|o  « 

at 

— 

•_ . 

H 

>*  1 

0 

CO  O   CO 

o 

O 

CO 

> 

LU 

CO     1 

1— 

0 

# 

LU 

2 

o 

< 

< 

•     1 

0 

ON    II 

or 

on 

LJ 

X 

O 

O     I 
1 

< 

LL 

O 

_ 

—  *- 

mm 

r> 

— 

IO 

Uj 

1 

ILU 

►-  CO  *- 

Z 

IN 

< 

K- 

co    1 

CO 

O 

Z  O  cc 

*— 

»» 

"^ 

< 

*-  1 

J 

t 

< 

<u     1 

< 

o 

X 

O 

O 

1 

O 

-4 

!• 

h- 

< 

57 

t 

Ml 

LU 

z 

or  1 

OJ 

co 

1 

mi 

•  « 

i-t 

o^ 

t- 

n 

LU     | 

•■ 

11 

ro 

z 

LL 

•• 

co 

O 

mm 

Q 

<r 

a 

1     1 

K 

0 

D 

t- 

Z 

UJ 

z 

o 

a 

H 

CO    I 

•— 

h-  C 

u 

I 

•-• 

C 

Q 

t— . 

a 

z   1 

m 

•— 1 

►- 

1- 

Z 

cc 

.— ■ 

z 

CO 

>— 

d 

<  1 

l 

cc 

LU 

•— 1 

1 

K 

X) 

II 

•— 

1 

< 

1 

~H 

CO 

CO 

< 

a  a 

CO 

t-  1 

l< 

II 

1 

c  OC 

Ul 

X 

Ul 

00 

0  1 

LL 

IZ 

-J 

-1 

.«  LL 

1— 

*V 

z 

LU 

LU     ( 

O 

a. 

-J 

< 

O  O    •» 

u 

LU 

< 

— 

*— 

— 

a      1 

1  O  O 

< 

> 

Z      IX 

a 

ISI 

~ 

X 

CO 

a: 

OT      | 

z 

IO 

0 

LU 

LU  Z  — 

r> 

— 

-J 

mm 

a 

LU 

c  1 

»— " 

z 

—  -1 

o 

-J 

< 

Z 

0 

CO 

0  1 

1/1 

1— • 

CO  Ui 

LU 

< 

> 

•— 1 

( 

0 

a. 

y 

— 

UJ 

CO 

1 

1 

1 

or 

>-  a. 

O  K 

1 

1 

t  1 

<T 

-J 

Z 

Q. 

— 

1 

U) 

< 

►1  0 

a 

Z  4 

1  * 

« 

<* 

» 

#  # 

3- 

> 

LU 

UL    Z 

a  lu 

0  0 

♦    I 

UJ  LU 

<M  >»• 

O  IA 

h»  0> 

CO   -4 

*\  ITl 

-4  O 

-^ 

l/l   r-4 

0 

00  -1 

1 

~o  ift 

LU 

m  m 

O 

m  (m 

— 1 

o»  h. 

r- 

a>  >t 

>»■ 

r»  *  * 


LU 

a 
a 


0000000000  00  000  000 
000000000000  000000 
-^rgro^ino^aocO-H icgro^irtot^oo 

000000000"-1"*•""•     m+   ~+   r+    r* I    -M    —< 

000000000000  000  000 


o  o 

*  I 

Ui  LU 

<\J   CO 

o  o 

f-  •»• 

00  ift 
rO  0> 

-J  -1   0*   -41 

o  in  >*■  o 

I  ooim   1 

LU   -O    >0    LU 

ao  ia  a>  o 
-o  c>  -o  o 
ft  o>  «*  «»• 

(Ni  ^  »-l  ^ 

o  r»  >r  <* 
ao    •    1    t 
O  f\l  co  as 
ut    I 
O 

o 
r» 
>*•  —  —  — 

•  or 

CO  o 

or 

or 

LU   •• 

"UJQ 
UJ  >  Z 
O  —  3 
-1  »-  O 
<  <   CO 

>  _l 

UJ    UJ    IO 

OKUQ 

uu      z.  z. 

t-   Z  IU   3 

0000 
a.      »  ca 

O  O  Z  CO 

O  Z  O  ►- 

Z3  O  CO 

O    I  >■ 

co  »*  -J 

<   •  z 

x  o»  < 

O  0" 

-•      -J 

--.   »    LU   > 

aj  u>  —1  ac 

I  C\j  _J  UJ 

<     v0    CL    (— 

m  m  m> 


3 
O 


<1J 
u 

S- 
Z3 
O 
CO 


or: 

I- 
or: 
o 
sz 
co 

CD 
Q. 

E 
(O 
X 


a. 


0 
z' 


33 


Example  0:  Computation  of  exp(8.0) 


This  program  involves  the  use  of  addition,  multiplication,  and 
division  by  computing  exp(8)  with  the  usual  power  series,  truncated 

o 

at  40  terms.  To  ten  digits  e  =  2980.957987,  and  at  the  end  statistical 
analysis  claims 

e8  e  (2980.9549...,  2980.961...) 


while  Interval  Analysis  gets 


e8  e  (2980.955   ,  2980.958) 


-  a  slightly  tighter  bound 
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Example  1:  Computation  of  cos(4tt) 

This  problem  involves  all  four  arithmetic  operations  (using 
again  the  usual  power  series)  and  is  not  well-disposed  to  statistical 
analysis  since  there  is  eventually  tremendous  cancellation  between 
terms.  Note  that  to  the  first  ten  terms  the  computed  sum  is  roughly 
-309,  and  the  eleventh  term  is  around  +387.6,  resulting  in  a  factor 
of  16  increase  in  variance.  Statistical  analysis  can  no  longer  give 
decent  results:  it  provides  the  bounds 

cos(4tt)  e  (-.4183606...,  2.4383515...) 

whereas  Interval  Analysis  gets 

cos(4tt)  e  (.7754054   ,  1.227057). 

Neither  of  these  results  are  good  but  the  statistical  bound  is  terrible. 
The  relative  error  bound  is  so  big  that  the  assumption  that  second  order 
effects  are  negligible  is  in  question. 
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Final  example:  Lanczos'  Algorithm 

Lanczos'  Algorithm  is  a  method  for  reducing  a  symmetric 
matrix  to  tridiagonal  form.  Without  re-orthogonal ization,  the  vectors 
of  the  orthogonal  transformation  which  it  generates  deviate  gradually 
from  orthonormality.  The  SNORTRAN  coding  of  the  algorithm  was  written 
before  the  problem  with  computed  zeroes  in  statistical  relative  error 
analysis  was  fully  appreciated,  otherwise  it  would  never  have  been  done. 
It  is  included  here  as  an  example  of  how  the  SNORTRAN  system  was  intended 
to  be  used  (the  PL/ I  preprocessor-output  is  included,  since  the 
algorithm  uses  so  many  SNORTRAN  features),  and  also  as  an  indicator  of 
how  severe  a  problem  the  irrepresentabil ity  of  zero  is.  The  program 
is  initially  given  v..  =(1,0,0,0,...)  and  is  asked  to  generate  63  other 
mutually  orthonormal  vectors.  Statistical  bounds  for  <v-,  ,v?>  are  huge 
because  the  result  is  so  close  to  zero,  and  during  computation  of 
<Vp,v3>  the  program  blows  up  because  a  zero  is  actually  encountered. 
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C0100 
00200 
00300 
00*00 
00500 
00600 
007OO 
00800 
00900 
01000 
01100 
01200 
01JOO 
01*00 
01500 
016OO 
017OO 
0180C 
01900 
C2000 
02100 

02300 
02400 
02500 
02600 
02700 
02800 
02900 
03000 
03100 
03200 
03300 
03400 
03500 
03600 
03700 
03800 
039OO 
04000 
04100 
04200 
04300 
0*400 
04500 
C4600 
04700 
04800 
C4900 
05000 
05100 
05200 
05300 
05400 
05500 
05o00 
05700 
05dOO 
05  900 
06000 


LANCZOS:       PRGC    OPT ICNS( MA  IN)     REORCEP; 

OECLAKE 

(  1  tJtK.I I,  IK, IX, IL,NLIMIT,N,P,NP)  FIXED  BIN, 
T   FLCAT  BINARYI53); 

SMEAR   0(64),  V(8,64),  rf(64),  ALPHA,  BETA,  GAMMA; 

ShEAR   VSAVE(3,64>,  IPVALlb),  TEMPU64),  TEMP2164); 

DECLARE 

STRING  (8)  ChAR(13)  I  N  IT  (  •  <V  I J  )  ,  V  (  J  )  >  •  ,  •  <V  (  J  )  ,  VU-  1 )  >•  , 
•<V<  J)  ,V(  J-2I>%  «<Vi  J)  ,VIJ-3I>'  ,,<V(J),V(J-4)>«, 
•<V(J),V(J-5)>»,«<V(J),V(J-6)>',,<V(J),V(J-7)>«); 

INITIALIZE  ROUNDED; 

N  =  8;  P=8; 
NP  =  N*P; 


/* 


INI T  V(l)  »/ 

CALL  SETC0NST4NT(V(1,1),SJRT( I . 0 000000000 COOOOOOOOCOOE 08 /NP J ) ; 
DC  <=1  TC  NP; 
VSAVE(1,K)  ,  VU,K)  =  Vll.l); 

END; 


/* 


/* 


Computation  of 
the  square  root 
of  the  smeared 
number  GAMMA 


INIT    Ull)    =   A   v(  1) 
DO     11  =  1     TO    NP;  TEMP1UI)     =     V(l,II); 

CALL    AMbLTtO, TEMPI); 


MAIN    LOOP 
DO    IX=2    TO    NP; 
IK    =    MODI lX-2,8)    +    l; 
DO    I  1=1    TO    NP;         TEMPKI  I) 
CALL    INNERIALPHA, TEMPI, U) ; 

DO    K=l    TO    NP; 

EVAL    W(K)    =    U(K)    - 

END; 
CALL     INNfJP<G&.MMA,  w,w)  ; 


END; 


*/ 


*/ 


=  V<  IK, II  ); 


END; 


ALPHA  *  V( IK, K) ; 


GAK.-1A.  Cof  PUTcD_V  A  LUl£  =  SCkT(  GAMMA.  CCMPGTED_  VALUE)  ; 
GAMMA. VARIANCE  =  GAMMA. VAR IANCE/2 ; 
T  =  GAMMA.IA_Lb;  GAMMA.  I  A_LB=RCUNDUCbN  I  SQRT  (  T)  ); 
T  =  GAMMA.I  A_Ob;  GAMMA.  I  A_UB=RGUNOUP   (SQRT(T)); 


/* 

/* 


1L  =  M0D(IK,8)  +  1; 

DO  K=l  TO  NP; 

EVAL  VUL.K)  =  W(K)  /  GAMMA; 

END; 

BETA  =  GAMMA; 

DO  11  =  1  TO  NP;     TEMPHII)  =  VtlL.II);    END; 
CALL  AMULTIU, TEMPI) ; 

DO  K=l  JD    NP; 

EVAL  U(K)  =  U(K)  -  BETA  *  V(IK.K); 

END; 

END  bF  MAIN  LOOP  CCKPUTATION  */ 

NO'*  CHECK  INNER  PRODUCTS...  */ 

IF  IX<8  1H2N  NLIM1T=IX;  ELSE  NLIMIT=b; 

PUT  SKIP  EOIT  I'J  =  ',IX)  (SKlP(2),A,F(3) )  ; 

DC  1=1  TO  NLI^IT; 

J  =  MCD(  I  L+ct-  I  ,  d  )  ♦  l; 

CO  11  =  1  TO  NP;  TEMPKII1  =  V(IL.II);   END; 

CU  11=1  TO  NP;  TEMP2III)  =  V(J,II);    END; 
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06100  CALL  INNEM1PVAL(  I  ) , TEMPI ,TEMP2 )  ; 

06200  CALL  SME ARED.PRI NT (STR ING ( I ) , IP VAL ( I  J ) ; 

06300  END; 

C6400  END; 

06500 

06600  AMOLT:   PK0C(U,V); 

06700  SMEAR  U(*l  »  V(*) ; 

06800  DECLARE  U,J,K)  FIXED  BIN; 

06900  K  =  0; 

07000  DO  J=l  TO  P; 

07100  DO  1=1  TO  N; 

07200  K  =  K  +  l; 

07300  EVAL   U(K)  =  4  *  V(KI; 

07^00  If  I>1  THEN   EVAL    U(KJ  =  U(KJ  -  V(K-l); 

07500  IF  KN  THEN   EVAL    U(KJ  =  UIK)  -  V(K+1); 

07600  IF  J>1  THEN   EVAL    U(K)  =  U(KJ  -  V(K-N); 

07700  IF  J<P  THEN   EVAL    UIK)  =  UIK)  -  V(K+N); 

07800  END; 

07900  END; 

08000  END  AMULT; 

08100 

C8200 

08300  INNER:   PROC ( I P, U, V) ; 

08400  SMEAR  IPt  U(*),  V(*); 

08500  DECLARE  K  FIXED  BIN; 

08600  IP  =  o; 

08700  DO  K=l  TO  NP ; 

08800  EVAL  IP  =  IP  ♦  U(K)  *  V(KJ; 

08900  END; 

09000  END  INNER; 

09100 

09200  END    LAivlCZOS; 


Final  Example:  Lanczos'  Algorithm  SNORTRAN  Source 
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MT  LEV  NT 


Preprocessor  Output  of  LANCZOS  Program 


1     o 

2    1   0 


3    1   0 


*    1   0 


5    1   0 


6    10 


0 
0 
0 
0 
0 
0 
0 

1 

0 
0 

0 


LANCZOS:   PROC 
DECLARE 


OPTIONSIMAIN)  REORDER; 

II,J.K,II,IK,IX,IL,NLIMIT,N,P,NP) 
T   FLOAT  BINARY(53); 


FIXED  BIN, 


DECLARE 

1  U(64)  LIKE  SMEARED.NUMBER, 
1  Vf8,6*)  LIKE  SMEARED_NUMBER, 
I  W(6*>  LIKE  SMEARED_NUMBER, 
1  ALPHA  LIKE  SMEARED_NUMBER , 
I  BETA  LIKE  SME ARED_NUMBER , 
1  GAMMA  LIKE  SMEARED_NUMBER; 


DECLARE 

i    VSAVE(8,64)    LIKE    SMEARED_NUMBER, 
I    IPVALI8)    LIKE    SMEARED_NUMBER, 
I    TEMPK64)    LIKE    SME ARED.NUMBER, 
I    TEMP2I64)    LIKE    SME ARED_NUMBER; 
DECLARE 

STRING    (8)    CHARI13)     IN  I T( • <V( J ) ,V{ J )> • , • <V  I  J I  ,  V  (  J-l  )>  •  , 
•<V{ J),V(J-2)>« f'<V(J) tVt J-3)>* ,»<VCJ),V( J-4)>' , 
,<V{J),V(J-5)>,,«<V(J),VlJ-6)>« ,,<V(J),V(J-7)>« ); 


DECLARE 
1 
2 
2 
2 


DECLARE 


SMEARED_NUMBER, 

COMPUTED_VALUE  FLOAT 

VARIANCE  FLOAT 

INTDATA, 

3    IA_LB  FLOAT 

3    IA_UB  FLOAT 


BINARYI21 ), 
BINARY<21), 

BINARYI21 ), 
BINARY121); 


/♦UNITS  BETA**(-2T)*/ 


/* 
/* 


INTERVAL 
ANALYSIS 


DATA 


*/ 
*/ 


I  REG(IO)  LIKE  SMEARED.NUMPER, 

8ETA_MT   FLOAT  BINARYI21)  INIT ( 1.0E-24B1 , 

HUGE      FLOAT  BINARYI21)  INIT < 1 .0E+88B ) ; 


/*l6**(-6)*/ 


SMEARED_ADD:   PRQC  (A,B,CI; 

/*   PERFORMS  ADDITION   A  =  B  +  C   FOR  SMEARED  NUMBERS   */ 
DECLARE  I  (A,B,CJ  LIKE  SMEARED.NUMBER , 
IT,T1,T2,T3)  FLOAT  BINARY153); 
T  =  B.COMPUTED_VALUE; 

A.COMPUTED_VALUE  =  ROUNDIT  ♦  C .COMPUTED_VALUE 1 ; 
T2  =  MAXSO(B. INTDATA); 
T3  =  MAXSQIC. INTDATA); 

CALL  INTERVALU.  INTDATA,  1  /*ADD*/,  B.  INTOATA,  C.  I  NTD  ATA)  ; 
Tl  *  MINSOU. INTDATA); 
IF  KT2+T3I/T1)  >  HUGE   THEN 

DO;   A. VARIANCE  *  HUGE;  RETURN;   ENO; 
A. VARIANCE  =  ROUNDUP!  T2/T1*8.VAR  I ANCE  ♦  T3/T l*C .VARI ANCE 
CALL  ADD_ROUNDING_ERROR(A); 
END  SMEARED_ADD; 


) ; 
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23 
24 

25 

26 
27 
28 
29 
30 

31 

32 

33 

34 
35 
36 
37 
38 
39 
40 
43 
44 
45 

46 

47 

48 
49 
50 
51 
52 
53 
54 

55 

56 

57 
58 
59 
60 
61 
62 
63 

64 

65 

66 
67 
68 
69 
70 

71 


1  0 

2  0 


1  0 

2  0 


1  0 

2  0 


SMEARED_NEGATE:  PROC  (A, B) ; 

/*   PERFORMS  NEGATION   A  »  -B   FOR  SMEARED  NUMBERS   */ 
DECLARE   1  (AtB)  LIKE  SMEAR ED.NUMBER, 

T  FLOAT  BINARY(21); 
A.COMPUTED_VALUE  =  -B.COMPUTED.VALUE; 
A. VARIANCE  =  B. VARIANCE; 
T  =  B.IA_LB; 
A.IA.LB  =  -B.IA_UB; 
A.IA_UB  =  -T; 
END  SMEARED.NEGATE; 

SMEARED.SUB:   PROC  (A,B,C); 

/*   PERFORMS  SUBTRACTION   A  =  B  -  C   FOR  SMEARED  NUMBERS   */ 
DECLARE  1  (A,B,C1  LIKE  SMEARED.NUMBER , 
<T,Tl,T2,T31  FLOAT  BINARYI53); 
T  =  B.COMPUTED_VALUE; 

A.COMPUTEO_VALUE  =  POUNDIT  -  C.COMPUTED_VALUE ) ; 
T2  =  MAXSQ(B.INTDATA) ; 
T3  =  MAXSQIC.INTDATAl; 

CALL  INTERVAL(A.INTDATA,2  /*SUB*/t  B. INTDATA, C. I NTOAT& )  ; 
Tl  =  MINSQ(A.INTDATA) ; 
IF  (|T2*T3»/TU  >  HUGE   THEN 

DO;   A. VARIANCE  =  HUGE;  RETURN;   END; 
A. VARIANCE  =  ROUNDUP!  T2/T1*B. VAR I ANCE  ♦  T3/T 1*C . VARI ANCE  ); 
CALL  ADD_RnuNDING_ERROR(A); 
END  SMEARED.SUB; 

SMEARED_MULT:   PROC(A,B,C); 

/*   PERFORMS  MULTIPLICATION   A  =  B  *  C   FOR  SMEARED  NUMBERS   */ 
DECLARE   1   (A,B,C>  LIKE  SMEARED_NUMBER t 

T  FLOAT  BINARYI53) ; 
T  =  B.COMPUTED_VALUE; 

A.COMPUTED_VALUE  =  ROUNDIT  *  C .COMPUTED_VALUE I ; 
CALL  INTERVAL(A.INTDATAt3  /*MULT*/ , B. I NTDATA, C. I NTD ATA ) ; 
T  =  B. VARIANCE; 

A. VARIANCE  =  ROUNDUPIT  +  C. VARIANCE); 
CALL  ADD_ROUNDING_ERR0R(AI; 
ENO  SMEARED.MULT; 

SMEARED.DIV:   PROCIA.B.C); 

/*   PERFORMS  DIVISION   A  *  B  /  C   FOR  SMEARED  NUMBERS   */ 
DECLARF   I   (A,B,C)  LIKE  SMEARED.NUMBER , 

T  FLOAT  BINARYI53I; 

T  «  B.COMPUTED_VALUE; 

A.COMPUTED_VALUE  =  ROUNDIT  /  C . COMPUTED_VALUE ) ; 

CALL  INTERVAL(A.INTDATA,4  /*DIV*/, B . INTDATA ,C . INTDATA) ; 

T  =  B. VARIANCE; 

A. VARIANCE    =    ROUNDUPIT    ♦    C. VARIANCE); 

CALL    ADD_ROUN0ING_ERROR(A) ; 

END    SMEARED_DIV. 


SETCONSTANT:     PROC(A,C); 

/*       SETS    SMEARED    NUMBER       A       TO    THE 
DECLARE        I       A       LIKE    SMEAR ED_NUMBER» 

C    FLOAT     BINARYC53); 
A.COMPUTED_VALUE    =    ROUND(C); 
A.IA_LB    =    ROUNDDOWN(C) ; 
A.IA_UB    =    POUNDUP(C); 
A. VARIANCE    =    OEOB; 
IF    C    -=    PRECISIONS, 21)     THEN 

CALL    ADD_ROUNDING_ERROP (A) ; 
END    SETCONSTANT; 


INITIAL    VALUE       C      */ 
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i  o 

2   0 


86  2   0 

87  2   0 

88  2   0 


2  0 

2  0 

1  0 

2  0 


INTERVAL:   PROC ( A, OP, B, C  ) ; 

/*   PERFORMS  INTERVAL  ARITHMETIC    A  =  B  OP  C    */ 
DECLARE  1  INTERVAL, 
2  IA.LB 
2  IA_UB 
1  (A,B,C)  LIKE  INTERVAL, 
(SL,S2,S3,S4)  FLOAT  BINARY(53), 
I0P(4)  LABEL, 
OP      FIXED  BINARYU5); 
SI  =  B.IA_LB;  S2  =  B.IA_UB;  S3  =  C.IA_LB;  S4  =  C.IA.UB; 
GO  TO  IOP(OP); 


MAXSQ 


FLOAT  BINARY(21), 
FLOAT  BINARY(21), 


IOP(l): 
A.IA_LB 
A.IA_UB 
RETURN; 

IOPC2): 
A.IA_LB 
A.IA_UB 
RETURN; 

IOP(3): 

A.IA.LB 


A.IA_UB  =  MAX( 


/*  ADDITION  */ 
ROUNDOOWNfSl  +  S3); 
ROUNDUP   (S2  ♦  S4); 


/*   SUBTRACTION   */ 
ROUNDDOWNIS1  -  S4); 
ROUNDUP   (S2  -  S3); 


/*  MULTIPLICATION  */ 
=  MIN(  ROUNDDOWN(Sl 
ROUNDDOWN( SI 
ROUNDDOWN(S2 
ROUNDDOWN(S2 
ROUNDUP  (SI 
ROUNDUP  (SI 
ROUNDUP  (S2 
ROUNDUP   (S2 


* 

S3), 

* 

S4)  , 

* 

S3)  , 

* 

S4) 

* 

S3), 

* 

S4), 

* 

S3), 

* 

S4) 

RETURN; 

I0PI4I: 

A.IA.LB 


A.IA_UB  *  MAX( 


/ 

S3), 

/ 

S4), 

/ 

S3), 

/ 

S4) 

/ 

S3), 

/ 

S4), 

/ 

S3), 

/ 

S4) 

); 


/*  DIVISION  */ 
MIN(  ROUNDDOWNtSl 
ROUNDDOWN(Sl 
ROUNDDOWN(S2 
ROUNDDOWN(S2 
ROUNDUP  (SI 
ROUNDUP  (SI 
ROUNDUP  (S2 
ROUNDUP  <S2 
RETURN; 

END    INTERVAL; 

PROC(A)       RETURNS(FLOAT    BINARY(53)); 
/*       RETURNS    MAX    VALUE    OF    A**2    OVER    INTERVAL    A 
DECLARE    1    A,       /*    INTERVAL       */ 

2    LB    FLOAT    BINARY(2l), 

2    UB    FLOAT    BINARY(2l), 

(T1.T2)    FLOAT    BINARY(53); 
Tl    ■    PRECISION(A.LB,53)**2; 
T2    =    PRECISION( A.U3,53)**2; 
Tl    =    MAX(T1,T2) ; 
RETURN(Tl); 
END    MAXSQ; 


*/ 
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99 
100 


101 
102 
105 
106 
107 
108 
109 


110 
HI 


1  0 

2  0 


1*0 


112 

2 

0 

113 

2 

0 

114 

2 

0 

115 

2 

0 

116 

2 

0 

117 

1 

0 

lie 

2 

0 

119 

2 

0 

120 

2 

0 

121 

2 

1 

122 

2 

1 

123 

2 

0 

12* 

2 

1 

125 

2 

1 

126 

2 

0 

127 

2 

0 

128 

1 

0 

129 

2 

0 

130 

2 

0 

131 

2 

0 

132 

2 

0 

133 

2 

0 

13* 

2 

0 

135 

2 

0 

136 

2 

1 

(MINSQt 

I 

I 

f 

t 


PROC(A)   RETURNSIFLOAT  BINARY153)); 
/*   RETURNS  MIN  VALUE  OF  A**2  OVER  INTERVAL  A 
DECLARE  1  A,   /*  INTERVAL   */ 
2  LB  FLOAT  BINARYI21I, 
2  UB  FLOAT  BINARY(21> , 

(Tl,T2)  FLOAT  BINARYI53); 


*/ 


IF 


/*       INTERVAL    CONTAINS    ZERO       */ 
RETURN(Tl);     END; 


<A.LB*A.UB    <=    OEOBl    THEN 
DO;       Ti    =    1E0B/HJGE; 
Tl    =    PRECISIONS. LB, 53)**2; 
T2    =    PRECISIONS. UB,53)**2; 
Tl    =    MIN(Tl,T2l; 
RETURNIT1 ); 
END    MINSQ; 


ADD_ROUNOING_ERR0R:     PROC(A); 

DECLARE    1    A    LIKE    SMEARED_NUMBER , 

(D,S)    FLOAT    BINARYI53); 
/*      COMPUTE  D   =    MIN(MANTISSASl)  */ 

/*  S    «       1    /    !2DI**2      /    3  */ 

/*  *    COMPUTATION    ROUNDOFF    ERROR    VARIANCE,    */ 

/*  UNITS    BETA»*(-2T).  */ 

IF    A. VARIANCE    >=    HUGE    I     A.COMPUTED_VALUE=0E0B    THEN    RETURN; 

D    =    MIN(MANTISSA(A.I A_LB) ,MANTISSA(A.IA_UB) ); 

S    =    IE0B    /    <D*0*12E0); 

A. VARIANCE  =  ROUNDUP ( A. VARI ANCE  ♦  01; 

END  add_rounding_error; 


MANTISSA:   PR0CIXI   RETURNS ( FLOAT  BINARYI53II; 
DECLARE  (X,Y)  FLOAT  BINARYI21); 
Y  =  ABS(X); 

DO  WHILE  (Y  >  1.0000000000000000000000E08); 

Y  *  Y  *  1.0000000000000000000000E-4B; 

end; 

do  while  (y  <  1.0000000000000000000000e-a8); 

y  =  y  *  l.0oooo0o0o0oooo00ooo0ooe+*8; 

END; 
RETURN(Y) ; 
END  MANTISSA; 

SMEARED_PRINT:   PROC(AA,A); 

DECLARE  I  A  LIKE  SMEARED_NUMBER , 
AA   CHAR!*), 

(S,S1*S2, ABS1,ABS2,ABS3,ABS*)  FLOAT  BINI53); 
S  =  SORTS. VARIANCE); 

51  =  -3.625  *  S;    /*  99.71  CONFIDENCE  LIMITS   */ 

52  =  +3.625  *  S; 

A8S3  =  A.COMPUTED_VALUE  *  (1  ♦  S1*BETA_MT); 
ABS*  =  A.COMPUTED_VALUE  *  (1  ♦  S2*BETA_MT); 
IF  A.COMPUTED_VALUE  <  0   THEN 

DO;   ABS1=ABS3;  ABS3=ABS*;  ABS*=ABS1:   END; 

PUT  FILE(SYSPRINT)  SKIP  EDIT 

J AA, 'COMPUTED  VALUE:* , A.COMPUTED_VALUE } 

ISKIP,COL(7),A,X(8),A,E(25,15,16H 

(  »3.625«SIGMA  BOUNDS  ON  RELATIVE  ERROR: *, '(•, SI , •   ,',S2, 
•  I  *  BETA**C-T)'  ) 

( SKIP, COL( 121, A, COL (50) ,A,2  ( E ( 25 , I  5, 16  I , A )  I 
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41    2   0 


42  1   0 

43  2   0 


*4 

2 

0 

<»5 

2 

0 

46 

2 

1 

47 

2 

1 

48 

2 

1 

49 

2 

I 

50 

2 

0 

51 

2 

0 

52 

2 

0 

53 

1 

0 

54 

2 

0 

55 

2 

0 

56 

2 

0 

60 

2 

0 

61 

1 

0 

62 

2 

0 

63 

2 

0 

64 

2 

0 

68 

2 

0 

69 

0 

71 

0 

72 

0 

73 

0 

74 

1 

75 

1 

76 

0 

79 

0 

80 

o. 

81 

1 

82 

1 

85 

1 

86 

1 

IMMPLIED  99.7*-C0NFIDENCE  BOUNDS:1,  •(',  ABS3,  •  ,*, 
ABS4,  '  )•   ) 
(SKIP, COL (12), A, COL (50), A,  2  (E( 25, 15, 16), A)) 

(MNTEPVAL  ANALYSIS  BOUNDS:1,  •  (  •  , A  .1 A_LB , •   f«,  A.IA_UB,«  l«| 
(SKIP,  COL  (12),  A,  COL  (50),  A,  2  (  E  (  16,6  ,7  )  ,  X.<  9)  ,A)  )  ; 

END  SMEARED_PRINT; 


ROUND:  PROC(X)  RETURNS! FLOAT  BINARY(21II; 
DECLARE  (X,HALF)  FLOAT  BINARY(53), 


FLOAT  BINARYI21), 
BIT  (64), 

BIT(8)  DEFINED  E  POSITION(l), 
BIT(l)  DEFINED  E  POS I TIONl 33 ) , 
BIT(64)  INIT( 

• 0O0O00000OO0OOOO00OO000O0OOO00O01 OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO' B) , 
FEXPT     BIT(8>  DEFINED  F  POSITION(l); 
E  =  UNSPEC(X); 


RX 

E 

EEXPT 
RBIT 
F 


IF 


/* 


ROUND 


*/ 


RBIT   THEN   DO; 

FEXPT  =  EEXPT; 
UNSPEC(HALF)  =  F; 
RX  *  X  ♦  HALF; 

END; 
ELSE    RX  =  X;         /*   TRUNCATE   */ 
RETURN(RX); 

END  ROUND; 

ROUNDUP:  PROC(X)  RETURNS  (FLOAT  BINARY(2D); 
DECLARE  X  FLOAT  BINARY(53), 

RX  FLOAT  BINARY(21) ; 
IF  X  >  OEOB   THEN  RETURN( ROUND ( X )) ; 

ELSE  DO;  RX=X;  RETURN(RX);  END; 
END  ROUNDUP; 

ROUNDDOWN:  PROC(X)  RETURNS (FLOAT  BINAPY(21)); 
DECLARE  X  FLOAT  BINARY(53), 

RX  FLOAT  BINARY(21); 
IF  X  <  OEOB   THEN  RETURN( ROUND( X) ) ; 

ELSE  DO;  RX=X;  RETURN(RX);  END; 
END  ROUNDDOWN; 

N  =  8;  P=8; 
NP  =  N*P; 

INIT  V( 1)  */ 

CALL  SETCONSTANT(V(l,i) ,SQRT( 1. 0000000000000000000000E03 /NP ) ) ; 
DO  K=l  TO  NP; 

VSAVE(1,K)  ,  V(1,K)  =  V(l,l); 
END; 

INIT  U( 1)  =  A  V( 1) 
DO  11  =  1  to  NP;    TEMPKII)  =  V(1,II);    END; 
CALL  AMULT(U, TEMPI) ; 

MAIN  LOOP 
DO  IX=2  TO  NP*. 
IK  *  MOD(IX-2,8)  ♦  l; 

DO  11  =  1  TO  NP;    TEMPKII)  *  V(IK,II);    END; 
CALL  INNER(ALPHA, TEMPI, U); 
DO  K=l  TO  NP; 


/* 


/* 


/* 


*/ 


*/ 
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187 
188 
189 
190 
191 
192 
193 
195 
197 
198 

199 

200 

201 

2  02 

205 
206 

207 
208 
209 


210 
212 
213 
214 
215 
218 
221 
222 
223 
224 
225 

226 


227 
228 
229 
230 
231 

2 
2 
2 
2 
2 

0 
0 
0 

1 
z 

232 
233 
2  34 

2 
2 
2 

z 

2 
2 

235 

2 

2 

2  36 

2 

2 

237 

2 

2 

238 
239 

?4,0 

2 
2 
2 

2 
1 
0 

CALL 
CALL 


CALL 


SMEARED_MULT{ REG< 1 ) , ALPHA, V(  IK,K  )  I ; 
SMEARED_SUB<W(K)fU(K),REG( 1) ); 

END; 
call  inner(gamma,w,w); 
gamma.computed_value=sqrt(gamma.computed_value); 

gamma. variance  =  gamma . vari ance/ 2; 

t=gamma.ia_lb; 

t=gamma.ia_ub; 

IL  =  M0D(IK,8)  +  l; 

DO  K=l  TO  MP; 

SMEARED_DIV(V( IL , K ) ,W( K) .GAMMA) ; 
END*. 


GAMMA. I A_LB=ROUNDDOWN<  SQRT (t ) ) 
GAMMA. IA_UB=ROUNDUP   (SQRT(T)) 


BETA  ■  GAMMA; 

DO  I  1=1  TO  NP;     TEMPHII)  ■ 
CALL  AMULT(UtTEMPl); 
DO  K=l  TO  NP; 


v(lLvll); 


end; 


CALL 
CALL 

/* 
/* 


SMEARED_MULT(REG(U  ,BETA,V(IK,K)1; 
SMEARED_SUB(U(K) ,U(K),REG( 1) ); 

END; 

END  OF  MAIN  LOOP  COMPUTATION 

NOW  CHECK  INNER  PRODUCTS... 
IF  IX<8  THEN  NLIMIT=IX;  ELSE  NLIMIT=8; 
PUT  SKIP  EDIT  I'J  =  »,IX)  (SKIP(2)f A,F<3) ); 

DO  1=1  TO  NLIMIT; 

J  =  M0D(IL+8-It8)  ♦  i; 

DO    11  =  1    TO    NP;    TEMPHII)    =    V(IL.II); 

DO    I  1  =  1    TO    NP;    TEMP2UI)    =    V(J,II); 

CALL     INNERl  IPVALd  ),  TEMPI  ,TEMP2  )  ; 

CALL    SMEARED_PRINT(STRING(I), IPVAL(I)) 

END; 


*/ 


END; 
END; 


AMULT; 


END; 

PROC(U.V) 


DECLARE 

1    U(*)    LIKE    SMEARED_NUMBER, 
I    V(*)    LIKE    smeared_number; 

DECLARE    (I,JfK)    FIXED    BIN; 

k  =  o; 

DO  J=l  TO  P; 
DO  1=1  TO  N; 
K  •   K  ♦    li 

CALL    SETCONSTANT(REG( I) ,4); 

CALL    SMEARED_MULTIU(K),REG( 1) ,V(K)); 

IF    I>1    THEN 
CALL     SMEARED_SUBIU(K),U(K) ,V(K-1) ); 

IF    KN    THEN 
CALL    SMEARED_SUB(UCK)  ,U(K),V(K«-1)); 

IF    J>1    THEN 
CALL     SMEARED_SUBCU(K)  ,UCK)  ,V(K-N)  )  ; 

IF    J<P    THEN 
CALL     SMEARED_SU3(U(K),U(K),V(K>NI ); 

END; 
END; 
END    AMULT; 
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SNORTRAN  System  Listings 

1.  SNORTRAN  Preprocessor  (SNOBOL)     pp. 49-55 

2.  SNORTRAN  System  Routines  (PL/I )    pp. 56-60 

(Referenced  as  file  ' SNORT. RND'  in  line  09000  of  preprocessor) 
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00100 
C0200 
C0300 
00400 
00500 
00600 
00700 
00800 
00900 
01000 
01100 
01200 
01300 
OIhOO 
01500 
01600 
01700 
01800 
01900 
02000 
02100 
02200 
02300 
02400 
02500 
02600 
02700 
02800 
02900 
03000 
03100 
03200 
03300 
03400 
035C0 
03600 
03700 
03800 
03900 
04000 
04100 
C4200 
04300 
04400 
04500 
04600 
04700 
04800 
04900 
05000 
05100 
05200 
05300 
05400 
05500 
C5600 
05700 
05800 
05900 
06000 


»-- 

* 
* 
* 
* 
* — 


-EJECT 

* 
* 


A  POSTERIORI  ROUNDOFF  ERROR  LANGUAGE  PREPROCESSOR 


6TRACE  =  1000 

DEFINE('FIL..(  VECTOR, STRING  II  ,VAL«  ) 

DEFINED PARSE. EX  PR (II t J t LOOK  AHEAD. TOKEN, LOOK AHEAD. SYMBGL • I 

OEFINEl 'SEMANTICS! PR C0NUM1 •» 

DEFINE!  'GETCHAR!  )«  ) 

DEF1NE(  'GETNCNBLANKd  •} 

DEFINE( 'SCAN! !• I 

DEFINE(  'ADD( !• I 

DEFINE! 'REMOVE! ) ■ ) 

OEFINE!  'EKIT(STRING)'  I 

DEFINE! • LASTEM IT (TARGET) XtY* J 

DEFINE! 'ALLOCATE  D.REG  I  II  •  I 

DEFINE! 'FREE. REGI ST ERl  II' I 

DEFINE!  'CLEAKREGS(  II'  1 

KEYWORD  =  ('INITIALIZE'  I  'SMEAR'  |  »EVAL'  I  'PRINT'!  .  KEY 

NREGS  =  10 

REGISTER  =  ARRAYtNREGSI 

SEMICOLON  =  1 

VARIABLE"  =  2 

EQUALS  =  3 

PLUS  =  4 

MINUS  =  5 

STAR  =  6 

SLASH  =  7 

LPAREN  =  8 

RPAREN  =  9 

FACTOR  =  10 

TERM  =  11 

EXPR  =  12 

ASSMT  =  13 

ACCEPT  =  14 

COMMA  =  20 


SLR(l)  EXPRESSION  PARSER  TABLES 

DEPTH  =  50 

STATESTACK  =  ARRAY(DEPTHI 

SYM30LSTACK  =  ARRAY(DEPTH) 

TOKENSTACK  =  ARRAY(DEPTH) 

STORED. DEPTH  =  5 

STORED.  TOKE;M  =  A  RR  AY  t  STOR  EC.  Cb'PTH  I 

STORED. SYMBOL  =  ARRAY { STOP  ED. DEPTH  I 

NPRJDS  =  12;   NTRANS  =  76;   NSTATES  =  23 

3, 1, 2, 4, 5, 8, 10, 11, 12, 2, d, 10, 11,2, 8, 10 ,11, 2,' 


TS  =  '2,13 


1,4,5,1,4,5,6,7 


•4,5,8,  10,  11,12,1,4,5,0,7,9,  uttJiitti 

'9, 1, 4, 5, 6, 7, 9, 4, 5, 9, 2, 8, 10, 2, 8, IU, 2, 8, 10, 11,' 

•2, 8,  10 ,11, 1 , 4, 5, 0,7, 9,  1 , 4, 5, 6, 7, 9,' 

TA  =  ' 2,3,4,0,5,6,7,0,9, 10,11,5,8,9, 12,5,8,9,13, 5, • 

•6, 7, 3, 9, 10, 14, -5, -5, -5, 15, lu, -5, -2, 17, 10,-6,-6,-6,15, 16, • 
•-o, -7, -7, -7, 15, io,- 7, 17, 18, 19, 5, 8, 20, 5, 8, 21, 5, 8, 9, 22,' 
•5,8,9,23,-3,-3,-3,15,16,-3,-4,-4,-4,15,16,-4,' 

FT  =  '1,3,4,5,0,12,16,20,0,27,33,36,42,48,51,54,57,61,0,0,' 
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06100  ♦  •0,65/71,' 

06200  LT  =  '2,3,4,11,-12,15,19,26,-10,32,35,41,47,50,53,56,60,64,-11' 

06300  +  ' ,-3,-9, 70,76,  • 

06400  LS  =  '14, 13, 12, 12, 12, 12, 12, 11,11, 11, 10, 10,' 

G6500  RL  =  '2,3,3,3, 1,2, 2, 3,3,1,3,1,' 

06600  TRANS. SYMBCL  =  ARRAY ( NTRANS) ;   F  ILL( .TRANS. SYMBGL , TS 1 

06700  TRANS. ACTION  =  ARRAY  1 NTRANS )  ;   F I LL( .TRANS. ACT  ION ,T A ) 

06800  FIRST. TRANS  =  ARR AY( NSTATES ) ;   F I LL( .F IRST. TRANS , FT) 

06900  LAST. TRANS   =  ARRAY( NSTATES) ;   F ILL( .LAST .TRANS , LT ) 

07000  LHS. SYMBOL   =  ARRAY ( NPRODS J ;    F  ILL (  .LHS . SYMBOL , LS ) 

07100  RHS. LENGTH   =  ARRAY( NPRODS ) ;    F 1LL< .RHS .LENGTH, RL ) 

07200  ♦  :(REAO) 

07300  * 

074C0  FILL  1=1+1 

07500  STRING      SPAN! • 0123456789- • )     .    VAL    ','    =     :F(RETURN) 

07600  ITEMUVECTOR,I)    =    CONVERT  (  VAL  ,.  I  NTEGER)     :(FILL) 

07700  * 

07800  * 

07900  READ     CARO  =  INPUT  :F(END) 

08000  RESUME   CARO  dP  KEYWORD  :S(PREPARE) 

08100  OUTPUT  =  DIFFER(TRIMICARD) )   CARD        :(READ) 

C8200  * 

08300  PREPARE    CARD       LEMPI     .    SNARK    = 

08400  OUTPUT    =    0 IFFERl TR IM( SNARK ) )     •     •    SNARK       :($KEY) 

08500  * 

08600  *   INITIALIZE   —   CUPY  IN  SYSTEM  DECLARES  AND  SUBROUTINES 

08700  *  DETERMINE  ALSO  ROUNDING  MODE. 

08800  * 

08900  INITIALIZE  SCAM) 

09000  INITFILE    =    'SNORT. RND • 

09100  INITFILE    =    I  DENT ( SCAN (),' CHOPPED • )     'SNORT. CHP' 

09200  INPUTI.INIT, INITFILE) 

09300  INITCOPY       OUTPUT    =    INIT  :S(INITCOPY) 

09400  SCAN!) 

09500  OUTPUT  =  NE ( SCAN  TYPE , SEM ICCLON )  •   /*****  ERROR  IN  INITIALIZE' 

09600  ♦  •  STATEMENT  *****/•         : (RESUME) 

09700  * 

09800  *   SMEAR  —   GENERATE  DECLARATIONS  FCR  SMEARED  VARIABLES 

09900  * 

10000  SMEAR    SCAN() 

10100  OUTPUT  =  •   DECLARE' 

10200  SMEARLOOP     V  =  SCAN(  ) 

10300  EQ(SCANTYPE, VARIABLE)  :F(SMERRCR) 

10400  SCANl  ) 

10500  EQ(SCANTYPE, COMMA)  :F(SMEARFIN) 

10600  OUTPUT  =  '       L  '  V  '  LIKE  SMEAREU_NUMB  ER ,  •   MSMEARLOOP) 

10700  SfEARFIN   EQ ( SC ANTYPE , SEM I COLGN )  :F(SMERROR) 

10800  SMcNC    OUTPUT  =  •       1  •  V  '  LIKE  SMEAREC_NUM3ER;  '   KRESUME) 

109J0  SKbkRCR  OUTPUT  =  •   /*****  ERKCR  IN  SMEAR  LIST.  END  FORCED  *****/• 

1100U  :  (RESUME) 

11100  • 

11200  *   PRINT   —   FCRMATTE>  PRINT  CF  SMEARED  VARIABLES 

11300  * 

11400  PRINT  SCANO 

11500  P>  LIST       V    =     SCAM  ) 

11600  E<*(  SCANTYPt  .VARIABLE)  :F(PRERROR) 

11700  OUTPUT    =    »       CALL     S ME  A RED_PRI NT ( • "    V    "• , "    V    ");" 

llbOO  SCAN(  ) 

■  v-0  EfcHSCANTYPE, COMMA!  :S(PRLIST) 

120  E0( SC^NTYPE, SEMICOLON)  :S(RESUME) 
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12100  PRERROR  OUTPUT  =  •   /*****  ERROR  IN  PRINT  LIST.   END  FORCED  *****/• 

12200  MRESUME) 

12300  -EJECT 

12400  * 

12500  *   EVAL     —   PARSE  EXPRESSICNS  OF  SMEARED  NUMBERS  AND 

12600  *  GENERATE  APPRCPRIATE  SUBROUTINE  CALLS. 

12700  * 

12800  EVAL     SCANO 

12900  PARSE. EXPRO  :(RESUME) 

13000  * 

13100  SEMANTICS  :($(»PRUU*  PRODNUM)) 

13200  * * 

13300  *   l:  S   ==>   <ASSIGNMENT>  ;  * 

13400  * * 

13500  PRODI  :(RETURN) 

1370C  *   2:  <ASSIGNMENT>   ==>   <VAR>   =   <EXPR>  * 

13800  * * 

13900  PRCD2    LASTEMIT(TGKENSTACK<STACKPTR>) 

14000  CLEARREGS1)  :(RETURN) 

14100  * * 

14200  *   3:  <EXPR>   ==>   <tXPR>   ♦   <TERM>  * 

14300  *■ * 

14400  PKLD3    FREE . REGl S TEk ( TOKENSTACK<STACKPTR> J 

14500  FREE.REGISTER(TCKENSTACK<STACKPTR  ♦  2>) 

14600  AREG  =  ALLOCAT ED .REG( ) 

14700  EMIT!  ■   CALL  SMEARE D_ADD( '  AREG  •  ,'  TOKENSTACK<STACKPTR> 

14800  +  •»■  TCKENSTACkXSTACKPTR  +  2>  •);•  ) 

14900  TOKENSTACK<STACKPTR>  =  AREG  : (RETURN) 

15000  * * 

15100  *   4:  <EXPR>   ==>   <EXPR>   -   <TERM>  * 

15200  * * 

15300  PRCD4    FREE.REGISTERJ TOKENST ACK<STACKPTR> ) 

15400  FREE.REGISTERiTOKENST ACK<STACKPTR  +  2> ) 

15500  AREG  =  ALLCCAT ED .REG ( ) 

15600  EMIT!  ■   CALL  SMEARED.SUB ( •  AREG  •,'  TOKENST ACK<STACKPTR> 

15700  ♦  *t'  TOKENSTACK<STACKPTR  +  2>  •);'  ) 

15800  TOKENSTACK<STACKPTR>  =  AREG  :  (RETURN) 

15900  * * 

16000  *   5:  <EXPR>   ==>   <TERM>  * 

16100  * * 

16200  PR0C5  MRETURN) 

16400  *   t>:  <EXPR>   =  =  >   +   <TERM>  * 

16500  * * 

16600  PFGCo  TUKENSTACK<STACKPTR>    =    TOKENSTACK<STACKPTR    ♦    1>       : (RETURN) 

16800  *       7:  <EXPR>       ==>       -       <T£RM>  * 

17000  PRCD7         FkEE. REGISTER!  TfiKENST  ACK<STAC*PTR    ♦    1>) 

17100  AREG    =    ALLGCATED.R^Gt ) 

17200  EMIT(     •       CALL    SME ARED_NEGATE ( •       AREG       •,' 

17300  +  TCKENSTAC«STACKPTR     ♦    1>     •);•     ) 

17400  TuKENSTACK<STaCKPTR>     =    AKEG  : (RETURN) 

17500  * * 

176C0  *   b:  <TERM>   ==>   <TERh>   *   <FACTOR>  * 

17700  * * 

17800  Pt-008  FKEE.KEGISTFR(TGKENSTACK<STACKPTR>) 

17900  FFEE.RcGISTER( TUKENSTACK<STACKPTK    +    2> ) 

18000  AREG    =    ALLCCATED.REGt I 
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18100  EMIT!  »   CALL  SMEA RED.MULT ( •  AREG  •»•  TOKENSTACK<STACKPTR> 

18200  +  •»•  TCKENSTACK<STACKPTk  +  2>  ■)!'  I 

18300  TOKENSTACK<STACKPTR>  =  AREG  '.(RETURN) 

18400  * * 

18500  *   9:  <TERM>   ==>   <TERM>   /   <FACT0R>  * 

18700  PR0D9    FkEE.REGISTER(TUKENSTACK<STACKPTR>) 

18800  FREE.P.EGISTER{TOKENSTACK<STACKPTR  ♦  2>) 

18900  AREG  =  ALLOCATED . REG( ) 

19000  EMITi  •   CALL  SMEARED_DI V ( •  AREG  ','  TOKENST ACK<STACKPTR> 

19100  +  '  t'  TOKENSTACK<STACKPTR  +  2>  •);•  ) 

19200  TOKENSTACK<STACKPTR>  =  AREG  :  (RETURN) 

19300  * * 

19400  *       10:  <TERM>       ==>      <FACTCR>  *    ' 

19500  * * 

19600  PRCD10  : (RETURN) 

19800  *   11:  <FACTOR>   ==>   (   <EXPR>   )  * 

19900  * * 

20000  PRGOU   TCKENSTACK<STACKPTR>  =  TQK£NSTACK<ST ACKPTR  +  1>   :  (RETURN) 

20100  * * 

20200  *  <FACTUR>   ==>   <VAR>  * 

20300  * * 

20400  PK0D12   EQIVARINFO.O)  :S(RETURN) 

20500  CCNSTANT.VAK   AREG  =  ALLCC ATEO.REG ( ) 

20600  EMIT(  «   CALL  SETCONST ANT ( •   AREG   •,' 

20700  ♦  TOKENSTACK<STACKPTR>   •);'  ) 

20800  TOKENSTACK<STACKPTR>  =  AREG  KRETURN) 

20900  * * 

21000  ALLOCATED. REG    I  =  LT(I.NREGS)   I  +  1  :F(ALL.USED) 

21100  REG1STER<I>  =  I DENT( REGI STER<  I> )  'X'     : F( ALLOCATED  .REG ) 

21200  ALLOCATE   ALLOCAT EC. REG  =  'REGC  I  •)•  KRETURN) 

21300  ALL. USED   OUTPUT  =  •   /*****  ALL  REGISTERS  USED  *****/'  :{ALLOCATE) 

21400  * 

21500  FkEE. REGISTER     REGNAME   «REGC    BREAKC)')  .  I    :F(RETURN) 

21600  REGISTER<I>  =  :<RETURN) 

21700  * 

21800  CLEARREGS        I  =  LT(I,NREGS)  I  ♦  1  :F(RETURM 

21900  REGISTER<I>  =  : (CLEARREGS) 

22000  * 

22100  EMIT     OUTPUT  =  D I FFER ( SA VED  .STR  ING )  SAVED. STRING 

22Z00  SAVED. STRING  =  STRING  KRETURN) 

22300  * 

22400  LASTEMIT   SAVED. STRING   (BREAM'(')  LEN(l))  .  X   BREAKC,')  = 

22500  ♦  X      TARGET  :F(ASSMT) 

22600  OUTPUT  =  SAVED. STRING 

22700  SAVEC. STRING  =  : (RETURN) 

22800  ASSMT    OUTPUT  =  •   •  TARGET  «  =  •  TOKENST ACK<ST ACKPTR  +  2>  •;• 

22900  ♦  :(RETURN) 

23JC0  -tJECT 

23100  * * 

23200  * 

23300  *  SLR(l)  EXPRESSION  PARSER 

2340U  * 

2  3  500  * * 

21600  PARSE. LXPk  ST  AT  EST ACK<1>    =     1 

23700  SYMBCl STACK<1>    =    0 

23d00  CURRENT. STATE    =     1 

23SJO  STACKPTk    =    2 

2*000  NfEDPEAD    ■     «YtS' 
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24100 
24200 

24300 
24400 
24500 
24600 
24700 
24800 
24900 
25000 
25100 
25^00 
25300 
25400 
25500 
25600 
25700 
25800 
25900 
26000 
26100 
26200 
26300 
26400 
26500 
26600 
267C0 
26800 
26900 
27000 
27100 
27200 
2730C 
27400 
27500 
27600 
27700 
2780C 
27900 
28000 
28100 
28200 
28300 
28400 
28500 
28600 
28700 
28800 
28900 
29000 
29100 
29200 
29300 
29400 
29500 
29600 
29700 
29800 
29900 
30000 


* 
PARSE 


FIND 
FCUND 


* 
SHIFT 


I     =    FIRST. TRANS<CURRENT.STATE> 

J    =       LAST.TRANS<CURRENT.STATE> 
EU(I.U) 

NLEDREAO    =    1  DENT ( NEEDRt AD )    'YES* 

LUCKAhEAU. TOKEN    =    SCANO 

LOOKAHEAC. SYMBOL    =    SCANTYPE 

EQILCUK  AHEAD.  SYMBOL,  T RANS . SYMBOL <1 >) 

I    =    LT(1,J)     I     ♦     1 

J  =  TRANS. ACTIGN<I> 

EQ(JtO) 

LT[J,0) 


S(REDUCE.LRO) 
S(FIND) 


S(FOUND) 
S(FINO)F(SYNTAXERR) 

S(ACCEPT) 
S(REDUCE.LRl) 


REDUCE. 
REDUCE. 
Rfc'DUCc 


GCTO 
PUSHGOT 


STATESTACK<STACKPTR>  =  J 

SYMBOLSTACK<STACKPTR>  =  LOOKAHEAD. SYMBOL 

TOKLNSTACK<STACKPTR>  =  LOOKAHEAD. TCKEN 

CURRENT. STATE  =  J 

STACKPTR  =  LT(STACKPTR, DEPTH)  STACKPTR  ♦  1   :S(PARSE) 

OUTPUT  =  •    /*******  STACK  OVERFLOW  ****»*/»  : (PARSE) 

LRO       NEEUREAD  =  'YES1  UREDUCE) 

LR1       NEEDREAD  = 

PROD  =  -J 

STACKPTR  =  STACKPTR  -  RHS .  l.ENGTH<PROD> 

CURRENT. STATE  =  STATESTACK<ST ACKPTR  -  1> 

SYMBOLSTACK<STACKPTR>  -  LHS. SYMBOL<PROD> 

I  =  FIRST. TRANS<CURRENT.STATE> 

J  =   LAST.TRA,>iS<CURRENT.STATE> 

EQ(SYMdLLSTACK<STACKPTR>,TRANS.SYMBUL<I>)        :S(PUSHGOTO) 

I  =  LTU,J>  I  +  1  :S(GOTO)F(SYNTAXERR) 

0    J  =  TRANS. ACTION<I> 

STATESTACK<STACKPTR>  =  J 

CURRENT. STATE  =  J 

SEMANTICSIPROD) 

STACKPTR  =  STACKPTR  +  1  MPARSE) 


* 
ACCEPT 

SYNTAXERR  OUTPUT  = 

* 

* 

-EJECT 

*        SCAN  — 

* 
* 

SCAN 


/*****  SYNTAX  ERROR 


:  (RETURN) 
*****/t   :(FRETURN) 


LEXICAL  ANALYZER.      RETURNED  VALUE  IS  TCKEN  STRING. 
ALSO  RETURNS    SCANTYPE  —  NUMERICAL  TOKEN  CODE 

VARINFO  —  EXTRA  INFC  FOR  <VAR>  TOKENS 


SCAN 
NEXTC 

GETNL1 


HAR  = 
NBLANKl ) 


:($(  'CLASS*  CLASSNJ) 


* * 

*        CLASS  1  —  LETTERS,  BREAK  CHARACTERS:   BUILD  VARIABLES 

ClASSl   SCANTYPE  =  VARIABLE;   VAPINFC  =  0 
CLASS1A  ADK);   GETCHAR() 

GE(CLASSN,l)  LE(CLASSN,2)  :S(CLASS1A) 

CLASSlb  IDENT(NEXTCHAR, '  •)  :F(CLASS1C) 

GET.\Cr^BLANK(  ) 
CLASSIC  Ew(CLASSN,b)  :F(OUTI 

PARENCCUNT    =     1 
CLASSIC    AuOO  J       GETCHAR(  ) 
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30100 

PARENCOGNT  =  E0(CLASSN»8)  PAPENCOUNT  • 

►  1  :S(CLASS1D) 

30200 

PARENCOUNT  =  EQ(CLASSN,9)  PARENCOUNT  - 

-  1  :F(CLA3S1D) 

30300 

EUPARENCOUNT.O) 

:F(CLASS1D) 

30400 
30500 
30600 
30700 
30800 

AOD(  ) 

:(OUT) 

* 

CLASS  2  —  DIGITS:   8UILD  CONSTANTS 

CLASS2 

SCANTYPE  =  VARIABLE;   VAR1NFG  =  1 

30900 

CLASS2A 

APC(J;   GETNONBLANM ) 

31000 

VARINFO  =  VARINFO  +  1 

31100 

EQ(CLASSN,21 

:S(CLASS2A) 

31200 

VARINFO  =  ICENT(NEXTCHAR,' .• )  100 

;S(CLASS2A) 

31300 

IDENTINEXTCHAR, 'E' ) 

:F(GIT) 

31400 

AODO;   GETCHARU 

31500 

G£(CLASSN,2)  L5(CLASSN,4) 

:F(OUT) 

31600 

CLASS2B 

ADOO;   GETChAR{) 

31700 
31800 
31900 
32000 

EG(CLASSN,2) 

:S(CLASS2B)F(0UT) 

* 

CLASShS  3-10  —  TERMINALS     3+  4-  5/ 

6*  7=  8(  9)  10; 

32100 

CLASS3 

ACO(l;   GETCHARl);   SCANTYPE  =  PLUS 

32200 

OIFFEKIUNARY)  EQ(CLASSN»2) 

:S(CLASS2)F(0UT) 

32300 

CLASS4 

ALOU;   GETCHAR();   SCANTYPE  =  MINUS 

32400 

3IFFERIUNARY)  £0 (CLASSN, 2 ) 

:SICLASS2)F(CUT) 

32500 

CLASS5 

ADD();   GETCHAR();   SCANTYPE  =  SLASH 

32600 

EU(CLASSN,6) 

:F(OUT) 

32700 

CLASS5A 

REMQVEO;   GETNONBLANKl) 

32800 

CLASS5B 

EU(CLASSN,6) 

:F(CLASS5A) 

32900 

REMOVED;   GETNONBLANKO 

33000 

EQ<CLASSN,5) 

:F(CLASS5B) 

33100 

REMOVE!)  . 

:(SCAN) 

33200 

CLASS6 

AOD( ), 

;   SCANTYPE  =  STAR 

;10UT) 

33300 

CLASS7 

AL>D(  )  . 

•   SCANTYPE  =  EQUALS 

:(0UT) 

33400 

CLASS8 

AUDI ) ! 

,   SCANTYPE  =  LPAREN 

:(OUT) 

33500 

CLASS9 

ACDU  , 

SCANTYPE  =  RPAKEN 

:(OUT) 

33600 

CLASS10 

ALU!  ) 

•  SCANTYPE  =  SEMICOLON 

: (OUT) 

33700 

CLASS11 

ACOU  . 

SCANTYPE  =  COMMA 

:<OUT) 

33800 

CLASS12 

IMPOSSIBLE 

33900 

CLASSO 

REMUVE( ) 

34000 

OUTPUT  =  "   /******  EFROR,  UNKNOWN  CHARACTER  ■" 

34100 

♦ 

NEXTCHAR  "•  ENCOUNTERED;   CELETED. 

******/."   :(SCAN) 

34200 

* 

34300 

OCT 

UNARY  = 

34400 

UNARY  =  EQ(SCANTYPE.ECUALS)  •!• 

34500 

UNARY  =  EO(SCANTYPE.LPAREN)  ■!• 

34600 
34700 
34800 

♦ 

J(RETURN) 

GETCHAR 

CARU  LEN( 1  )  .  NEXTCHAR 

:F(MCRE) 

34900 

CLASSN  =  IOENTtNEXTCHAk, •  ')  12 

:S(RETURN) 

3500G 

CLASSN'  =  IDtNKNEXTCHAR,  •         •)  12 

:S(RtTURN) 

35100 

CLASSN  =  ICENTlNtXTUUk,  •  ,  •)  11 

:S(KETURN) 

35200 

CLASSN  =  REPLACE (NExTCHAk, 

35300 

♦ 

•AdCDEFGHIJK  LMNOPl.RSHlVWXYZa$_»0  12345676  9*-/*=  (  )  ;'  , 

35400 

♦ 

•OOOOOOUUOOOOOOJJCOOOCClICOOOOOOI  111111111^34  56  789') 

3550U 

CLASSN  =  INTLGER(CLASSN)  CLNVERT < CLASSN, .  INT EGE R )  ♦  I 

35600 

♦■ 

:S(RETURN) 

3570U 

CLASSN  =  0 

:  (RETURN) 

35800 

MIRE 

CARL  =  INPUT 

:S(GETCHAR) 

j  WOO 

OUTPUT  ■  •  /♦*«***  EOF  HIT  *******/• 

3o000 

CLASSf 

i    =  10 

'55 


36100 

NfcXTCHAR  =  ' ;• 

36200 

* 

36300 

&FTNCNBLANK      GETCHAR1 ) 

36*00 

EQJCLASSN,  12) 

36500 

CAPO  L£N(1)  = 

36600 

* 

36700 

ADO 

SCAN  =  SCAN  NEXTCHAR 

36800 

CARD  LEN(l)  = 

36900 

* 

37000 

REMOVE 

CARD  LENll)  = 

37100 

-EJECT 

37200 

END 

1RETURN) 


F(RETURN) 
(GETNONBLANK) 


(RETURN) 
IRETURN) 


V 
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00100 
00200 
00300 
00400 
00500 
00600 
00700 
00800 
00900 
01000 
01100 
01200 
01300 
01400 
01500 
01600 
01700 
01800 
01900 
02000 
02100 
OZZOO 
02300 
02400 
02500 
02600 
02700 
02800 
02900 
03000 
03100 
03200 
03300 
03400 
03500 
03600 
03700 
03800 
03900 
04000 
04100 
C42U0 
04300 
04400 
04500 
04600 
04700 
C4800 
04900 
05000 
05100 
05200 
C5J00 
05400 
05500 
05600 
05700 
05dOO 
05900 
06000 


DECLARE 
1 
2 
2 
2 


DECLARE 


SNORTRAN  System  Routines 


SMEARED_NUMBER, 

COMPUT£U_VALUE   FLCAT 

VARIANCE         FLCAT 

INTDATA, 

3    IA_LB        FLOAT 

3    IA_UB        FLOAT 


8INAkY<21) , 
BINARY(21), 

BINARY121) , 
BINARY(21 ); 


/♦UNITS  BETA**(-2T)*/ 


/*  INTERVAL 

/*   ANALYSIS  DATA 


*/ 
*/ 


1  REG(IO)  LIKE  SMEARED_NUMBER, 

BETA_MT   FLOAT  3INAkY<21)  IN  IT ( 1 .0E-24B ) , 

HUGE      FLOAT  BINARY(2i)  INI T < l.OE +888 ) ; 


/*16**(-6)»/ 


SMEAkfcO.ADD:        PRDC     (A,B,C); 

/*       PERFORMS    ADDITION       A    =    B    ♦    C       FOR    SMEARED    NUMBERS      */ 
DECLARE    1     (A,b,C)    LIKE    SME ARED_NUMeER, 
(T,T1,T2,T3)    FLOAT    BINARY(53); 
T    =    B.CCMPUTED.VALUE; 

A.CUMPOTED_VALUE  =  RCUND(T  +  C.CCMPUTED_VALUE ) ; 
T2  =  MAXSQ(B. INTDATA) ; 
T3  =  MAXSUIC. INTDATA) ; 

CALL  INTERVALU. INTDATA, 1  /*ADC*/,  B . INTDATA, C . I NTDATA ) ; 
Tl  =  MINSU(A. INTDATA) ; 
IF  UT2*T3)/T1)  >  HUGE   THEN 

DC;   A. VARIANCE  =  HUGE;  RETURN;   END; 
A. VARIANCE  =  ROUNDUP!  T2/T1*B. VAR I ANCE  ♦  T3/T1*C. VAR IANCE  ); 
CALL  AOO_RCUND1NG_EPRGR(A) ; 
END  SMEARED.ADD; 

SMEARED.NEGATE:    PROC    <A,B); 

/*   PERFORMS  NEGATION   A  =  -B   FOR  SMEARED  NUMBERS   */ 
DECLARE   1  (A,BI  LIKE  SME ARED.NUMSER , 

T  FLCAT  BINARYI2U; 
A.COMPUTED_VALUE  =  -b.COMPUTEC.VALUE ; 
A. VARIANCE  =  B. VARIANCE; 
T  =  B.IA_L3; 
A.IA_LB  =  -B.IA.UB, 
A. IA_UB  =  -T; 
END  SMEARED_NEGATE; 

SMEARED_SUB:   PRCC  (A,B,C»; 

/*   PERFORMS  SUBTRACTION   A  =  B  -  C   FOR  SMEARED  NUMBERS   */ 
DECLARE  1  (A,B,C)  LIKE  SME ARED_NUM8ER , 
(T,T1,T2,T3)  FLOAT  BINARY(53); 
T  =  B.CCMPUTED_VALUE; 

A.COMPUTED_VALUE  =  RCUNDIT  -  C.CGMPUTED.VALUE) ; 
T2  =  HAXSQIB. INTDATA) ; 
T3  =  MAXSy(C. INTDATA) ; 

CALL  INTEKVAUA.  INTDATA, 2  /*SUb*/»  b  .  INTDATA  ,C.  I  NTDATA  )  ; 
Tl  =  MI NSQI A. INTDATA) ; 
IF  ( (T2+T3J/TI )  >  HUGE   THEN 

DO;   A. VARIANCE  =  HUGE;  RETURN;   END; 
A. VARIANCE  ■  ROUNDUP!  T2/T  l*B . VAR I ANCE  ♦  T3/T 1 *C . VARI ANCE  H 
CALL  AUU.RCUNO ING_ERRLRIA) ; 
END  Sf'EARED_SUb: 


SME  AFtD_MULT  :   PfcCC  (  A  ,  B  ,C  )  ; 

/•   PLKFORfS  MULTIPLICATION   A  =  B  ♦  C   FOR 
ObCLAKL   1   (A,B,C)  L  IKt  SME AKED_NUMBER, 


SMEARED  NUMBERS   */ 
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06100 
06200 
0o300 
06400 
06500 
06600 
06700 
06800 
06900 
07000 
07100 
07200 
07300 
07400 
07500 
07600 
07700 
07800 
07900 
08000 
C8100 
08200 
08300 
08400 
08500 
06600 
08700 
C8800 
08900 
09000 
09100 
09200 
09300 
09400 
09500 
09600 
09700 
09800 
09900 
10000 
10100 
10200 
10300 
1040C 
10500 
10600 
10700 
10800 
10900 
11000 
11100 
11200 
11300 
11400 
11500 
11600 
11700 
11800 
11900 
12000 


T  FLOAT  BINARY153); 
T  =  B.COMPUTFD_VALUE; 

A.CJMPUTED.VALUE  =  RCUNDIT  *  C  .CCMPUTEO_VALUE )  ; 
CALL  INTERVAL<A.INTDATA,3  /*MULT*/ ,b. INTDATA , C . I NTGATA ) ; 
T  =  6. VARIANCE; 

A.VARIANCE--  RCUNDUPIT  ♦  C. VARIANCE); 
CALL  AOD_ROUNDING_ERRUR(A) ; 
END  SMEARED_MULT; 


SMEARED_DTV:   PROC(A,B,C); 

/*   PERFORMS  DIVISION   A  =  B  /  C   FOR  SMEARED  NUMBERS 
DECLARE   1   IA,B,C)  LIKE  SMEARED_NUMBER , 

T  FLCAT  61NARY(53); 
T  =  B.CUMPUTED_VALUE; 

A.COMPUTED.VALUE    =    RCUNDtT    /    C .COMPUTED_VALUE ) ; 
CALL     INTERVriL(A.INTDATA,4    /*C  1 V/ ,  B.  INTOATA,  C  INTDAT  A)  ; 
T    ■    B. VARIANCE; 

A. VARIANCE    =    RCUNDUP(T    ♦    C. VARIANCE); 
CALL    ADD_ROUNDING_ERR(JR(A)  ; 
E^D    SMEARED_DIV; 


»/ 


INITIAL  VALUE 


SETCCNSTANT:  PRGL(A,C); 

/*   SETS  SMEARED  NUMBER   A   TO  THE 
DECLARE   1   A   LIKE  S MEARED_NUMtiER , 

C  FLOAT  BINAPY(53); 
A.COMPUTED.VALUE  =  RCUNDIC); 
A.IA_LB  =  ROUNDDCnN(C); 
A.IA_Ud  =  RCUNDUP(C); 
A. VARIANCE  =  0E03; 
IF  C  -•=  PRECISIOMC21)  THEN 

CALL  ADD_ROUNDING_ERROR(A) ; 
END  SETCONSTANT; 


INTERVAL:   PROC ( A ,OP, B , C ) ; 

/*   PERFORMS  INTERVAL  ARITHMETIC    A  =  B  OP  C 
DECLARE  1  INTERVAL, 

2  IA  LB      FLCAT  BINARY(21) , 
2  IA_UB      FLCAT  BINARY<21>, 
I  (A,B,C)  LIKE  INTERVAL, 
(S1,S2,S3,S4)  FLOAT  BINAkY(53), 
I0P(4)  LABEL, 
OP      FIXED  BINARY(15); 
Si  =  B.IA_LB;  S2  =  B.IA_UB;  S3  =  c.ia_lb; 
GO  TO  ICP(CP) ; 


*/ 


*/ 


S4  =  c.ia.ub; 


IOP(  1): 
A.I A_Lb 
A.IA_UB 
RETURN; 

I0P(2»: 

A. IA_LB 
A.IA_UB 

return; 

IOP43)  : 
A.IA_LB 


/*  ADDITION  */ 
ROUNDDCWNlSl  ♦  S3)  ; 
ROUNDUP   (S2  *■    S4)  ; 


/*   SUBTRACTION   */ 
ROUNDDQWNISl  -  S4) ; 
ROUNDUP   (S2  -  S3) ; 


/*   MULTIPLlCATlbN   */ 
MINI    RGUNDLUWNtSi  *  S3), 
RCUNUU0WNCS1  *  S4 ) , 
R0UNCDGWN1S2  *  S3), 
R0UNLCJkN(S2  *  S4) 
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12100 
12200 
12300 
12400 

12500 
12600 
12700 
12800 
12900 
13000 
13100 
13200 
13300 
13400 
13500 
13600 
13700 
13300 
13900 
14000 
14100 
14200 
14300 
14400 
14500 
14600 
14700 
14800 
14900 
15000 
15100 
15200 
15300 
15400 
15500 
15600 
15700 
15800 
15900 
16000 
16100 
16200 
16300 
16400 
16500 
16600 
167C0 
16800 
16900 
17000 
17100 
17200 
17300 
17400 
17500 
17600 
17700 
17b00 
179  00 
16J00 


PAXSQ: 


MINSU: 


A.IA_U8  = 

MAX( 

ROUNDUP   (Si 

* 

S3), 

ROUNDUP   (SI 

* 

S4)  , 

RUUNDUP   (S2 

* 

S3), 

ROUNDUP   (S2 

* 

S4) 

RETURN; 

I0P(4): 

/* 

DIVISION   */ 

A.IA_LB  = 

MIM 

ROUNDDQwN(Sl 

/ 

S3), 

R0UNDD0WN(S1 

/ 

S4), 

RCUNDD0WN(S2 

/ 

S3), 

R0UNDD0WiM(S2 

/ 

S4) 

A.IA.UB  = 

WAX( 

ROUNDUP   (SI 

/ 

S3), 

ROUNDUP   (SI 

/ 

S4), 

ROUNDUP   (S2 

/ 

S3), 

ROUNDUP   (S2 

/ 

S4) 

RETURN; 

END  INTERVAL; 

PROC(A)   RETURNSIFLOAT  BINARY(53)); 
/*   RETURNS  MAX  VALUE  OF  A**2  OVER  INTERVAL  A   */ 
DECLARE  1  A,   /*  INTERVAL   */ 

2  LB  FLOAT  BINARYI21), 

2  UB  FLOAT  BINARY(21) , 

(T1,T2)  FLOAT  BINARY(53) ; 
Tl  =  PRECISICM A.LB,53)**2; 
T2  =  PRECISICN(A.UB,53)**2;  . 
Tl  =  NAX(Tl,T2); 
RETURN(Tl); 
END  MAXSQ; 

PROCU)   RETURNS(FLCAT  BINARY(53)); 
/*   RETURNS  MIN  VALUE  OF  A**2  OVER  INTERVAL  A   */ 
DECLARE  1  A,   /*  INTERVAL   */ 
2  LB  FLOAT  BINARY121) , 
2  UB  FLOAT  B1NARY(21), 

(T1,T2)  FLOAT  BINARY153) ; 


IF 


(A.LB*A.UB  <=  OEOB)  THEN 
DO;   Tl  =  IEOB/HUGE; 

=  PRECISION! A. LB, 53)**2; 

=  PRECISICN(A.UB,53)**2; 

=    MIN(T1,T2) ; 
RETUPM(Tl) ; 
END    MINSO; 


/*       INTERVAL    CONTAINS    ZERO      ♦/ 
RETURN(Tl);    END; 


Tl 
T2 

Tl 


iDD_RCJNMNG„ERRCR:    PROC(A); 

DECLARE     I     A    LIKE     SMEAKED_NUM6 ER , 

(D,S)     FLCAT    dINARY(53); 
/♦       COMPUTE  D    =    MIN(MANTISSAU))  »/ 

/*  S    =       1    /    (2D)**2       /    3  */ 

/*  =    COMPUTATION    RCUNOOFF    ERROR    VARIANCE,    */ 

/•  UNITS    BETA**(-2T).  */ 

IF     A. VARIANCE     >=     HUGE     |     A . COMPUT EO_VALUE=OEOB    THEN    RETURN; 

0    =    MINIMANT ISSA(A.IA_L8) , MANT I S SA ( A. I A_UB) ) ; 

S    =    1E03    /    (C*i)*12E0)  ; 

A. VARIANCE  =  kUUNQUPl  A.VAK  IANCE  *•  C); 

ENC  ADL_kOU>iDIN'-,_tKRCK; 
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18100 
18200 
18300 
18400 
18500 
18600 
1870C 
18800 
18900 
19000 
19100 
19200 
19300 
19V00 
19500 
19600 
19700 
19600 
19900 
20000 
2010C 
20200 
2030C 
20400 
20500 
20600 
20700 
20800 
20900 
21000 
21100 
21200 
21300 
21400 
21500 
21600 
21700 
21600 
21900 
22000 
22100 
22200 
22300 
22400 
22500 
22600 
22/00 
22800 
22900 
23000 
23100 
23200 
23300 
23400 
23500 
23600 
23700 
23d00 
23900 
24000 


MANTISSA:   PKOC(X)   RE  TURNS ( FLOAT  ETNARY(53)); 
DLCLARE  (X,Y)  FLOAT  bINARY(21); 
Y  =  AdS(X); 

00  WHILE  (Y  >  l.OGOOOOOOUOOOOOOOOOOOCOEOB) ; 

/  =  Y  *  1.0000000000000000000000E-4B; 

END; 

DO  WHILE  (Y  <  1.00000O0OCO0000UO00OO00E-4B) ; 

Y  =  Y  *  1.0000000COOCOOOOOOOOOOOE+4B; 

END; 
RETURN(Y); 
END  MANTISSA; 

SMEARED_PRINT:       PROC(AA.A); 

DECLARE    1    A    LIKE    SHE ARED.NUMBER, 

AA       CHAR(*), 

(S,S1,S2,ABS1,ABS2,ABS3,ABS4)  FLOAT  BIN(53); 
S  =  SQRTIA. VARIANCE); 

Si  =  -3.625  *  S;    /*  99.7%  CONFIDENCE  LIMITS   */ 
S2  =  +3.625  *  S; 

ABS3  =  A.COMPUTEO_VfiLUE  *  (1  +  S1*6ETA_MT); 
ABS4  =  A.CCMPUTEO_VALUE  *  (1  ♦  S2*bETA_MT); 
IF  A.COMPUTED_VALUE  <  0   THEN 

DO;   ABS1=ABS3;  ABS3=A8S4;  A6S4=ABS1;   END; 

PUT  FILE(SYSPRINT)  SKIP  EGIT 

(AA, 'COMPUTED  VALUE: • , A.COMPGTEO.VALUE ) 

(SKIP,C0L17),A,X(8),A,E(25.15,16)) 

{•3.625*SIGMA  BOUNDS  ON  RELATIVE  ERROR: «,•(', SI , •   ,',S2, 
•  )  *  BETA**(-T)«  ) 

(SKIP,CGL(12),A,C0L(50),A,2  (E{25,15,16),AJJ 


ABS4,  «  )•   ) 

(SKIP, COL  112), A, COL (50), A, 2  (E(25,15,16),A)I 

(•INTERVAL  ANALYSIS  BOUNDS:',  ■ ( • , A. I A_LB , •   ,',  A.IA_U8,'  )') 
(SKIP,CGL(12),A,CGL(5G),A,2  (E(16,6,7),X(9),A)); 


END  SMEAREC.PRINT; 


ROUND:  PROC(X)  RE  TURNS ( FLOAT  BINARY(21J); 
DECLARE  (X,HALF)  FLOAT  BINwRY(5i), 


FLOAT  8INARY(21), 
8IT(c>4) , 

3IT(8)  DEFINED  E  PUSITICN(l), 
3IT(1)  DEFINED  E  POS IT  I 0N( 33) , 
BIT(64)  IMT( 

•OOOOUOUOOOOOOOOOOOO 000000 000000 01 000000 00000000 00 0000 0000 30000 GO' B)  , 
FEXPT     BITlo)  DtrlNEC  F  POSITIGN(l); 
E  =  UNSPEC(X) ; 


RX 

E 

EEXPT 

RBIT 

F 


IF 


/* 


ROUND 


*/ 


RBIT   THEN   DO; 

FEXPT    =    EEXPT; 
UNSPEC(HALF)     =    F; 
RX    =    X    ♦    HALF; 
ENf  ; 
ELSE  KX    =    X  ;  /*       TRUNCATE       */ 
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24100  RETURN(RX); 

24200  ENO  RLUNO; 

24300 

24400     ROUNDUP:  PrtOC(X)  RETURNSC FLCAT  BINARY(21)1; 

24300  DECLARE  X  FLCAT  bINARY(53), 

24600  RX  FLCAT  BINARY(2ll; 

24700  IF  X  >  CEOB   THEN  RETURN ( RCUND (X )) ; 

24800  ELSE  00;  RX=X;  RETURN(RX);  END; 

24900  END  RCUNDUP; 

2500C 

25100  ROUNDDCJWN:    PROC(X)    RETURNS  ( FLOAT    EINARY(2D); 

25200  DECLARE    X    FLCAT    BINARY(53J, 

25300  RX  FLOAT  3INARY(21); 

25400  IF  X  <  OEOB   THEN  RETURN ( R0UNC( X )) ; 

25500  ELSE  DC;  RX=X;  FETURN(RX);  END; 

25600  END  RCUNODOWN; 
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7.  CONCLUSIONS 

The  important  conclusions  of  this  paper  are  first,  that 
accumulated  roundoff  errors  tend  to  cluster  towards  the  center  of  their 
bounds  like  normal  densities;  second,  that  high-confidence  bounds  can 
be  proved;  third,  that  automated  statistical  analysis  is  not  worthwhile. 

The  last  conclusion  is  not  obvious.  Naturally  the  statistical 
model  will  generate  better  bounds  than  the  worst-case  model  of  Wilkinson. 
But  beyond  that  no  conclusions  are  obvious.  Interval  Analysis  treats 
error  densities  as  uniform,  although  the  Central  Limit  Theorem  tells  us 
they  are  almost  normal.  We  would  expect  better  results  from  statistical 
analysis  because  it  not  only  combines  intervals  (read  "probability 
densities  on  intervals"),  it  also  retains  information  about  the  re- 
sultant error  distribution  on  the  new  interval.  But  statistical  analysis 
breaks  down  because  relative  errors  cannot  handle  cancellation  or 
numbers  near  zero,  and  correlation  is  almost  impossible  to  monitor 
practically. 

One  possible  exit  from  the  irrepresentability  of  zero  is  to 
drop  relative  errors  altogether  and  develop  a  statistical  theory  of 
absolute  errors.  Observations  2  and  3  of  Section  2  provide  most  of  the 
necessary  theory  if  we  now  treat  numbers  as  random  variables  instead  of 
the  number's  relative  errors.  Zero  then  becomes  representable.  How- 
ever, it  is  not  clear  that  an  analogue  of  Theorem  5  in  Section  4  can  be 
proved  for  absolute  analysis  since  multiplication  and  division  no  longer 
involve  simple  convolution,  and  besides,  the  problems  of  correlated 
errors  cannot  be  eliminated.  Statistical  error  analysis  in  general 
seems  of  questionable  use. 
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